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Permutohedral  aggregates  are  hierarchical  arrangements  of  cells.  Special  cases 
include  hexagonal  aggregates  and  truncated  octahedral  aggregates  in  dimensions  2 
and  3,  respectively.  These  cells  can  be  efficiently  addressed  using  a  scheme  known  as 
generalized  balanced  ternary,  or  GBT.  The  cells  of  a  Z)-dimensional  permutohedral 
aggregate  are  addressed,  or  indexed,  using  finite  strings  of  a  set  of  n  =  2D+l  - 1  digits. 
The  objects  of  interest  in  this  work  are  functions  whose  domains  are  permutohedral 
aggregates. 

The  generalized  matrix  product,  also  known  as  the  p-product,  is  a  generaliza- 
tion of  standard  matrix  multiplication. 

The  present  work  defines  the  discrete  Fourier  transform  and  the  p-product  for 
permutohedral  aggregates,  and  derives  a  radix-n,  decimation-in-space,  fast  Fourier 
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transform  (FFT)  for  permutohedral  aggregates  in  terms  of  the  p-product.  Data 
reordering  (also  known  as  shuffle  permutations)  is  generally  associated  with  FFT 
algorithms.  However,  use  of  the  p-product  makes  data  reordering  unnecessary. 
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CHAPTER  1 
INTRODUCTION 


The  goal  of  this  dissertation  is  to  describe  algorithms  that  efficiently  compute 
the  discrete  Fourier  transform  and  the  corresponding  inverse  transform  for  functions 
with  discrete  domains  known  as  permutohedral  aggregates.  The  Fourier  transform  is 
one  of  the  most  fundamental  transforms  in  signal  analysis  and  image  processing.  This 
transform  extracts  the  frequency  content  of  a  function.  Permutohedral  aggregates  are 
hierarchical  arrangements  of  certain  points  in  /^-dimensional  Euclidean  space.  These 
aggregates  are  of  interest  because  they  possess  several  useful  properties: 

1.  they  can  be  addressed  efficiently  using  finite  strings  of  digits; 

2.  they  have  a  ring  structure  and  addition  and  multiplication  can  be  done  entirely 
in  terms  of  bit  string  operations; 

3.  in  dimension  two,  aggregates  approximate  the  region  of  light-sensing  cells  in 
the  human  visual  system; 

4.  in  dimensions  two  and  three,  aggregates  represent  optimal  sampling  strategies 
in  the  spatial  domain. 

Algebraically,  A-level  permutohedral  aggregates  are  sets  of  coset  representa- 
tives of  quotient  lattices  of  the  form  L/NXL,  where  L  is  a  /^-dimensional  lattice 
having  a  ring  structure,  N  is  an  endomorphism  of  L,  A  is  a  positive  integer,  and 
NXL  =  NX(L).  These  quotient  lattices  consist  of  nx  elements,  where  n  =  2D+1  —  1. 


The  Voronoi  cell  of  a  point  x  in  a  lattice  /  is  the  set  {y  €  RD  :  \y  —  x\  < 
\y  —  2 1,  for  all  z  6  /}.  The  Voronoi  regions  or  cells  associated  with  permutohe- 
dral  aggregates  are  known  as  permutohedra.  From  an  image  processing  perspective, 
the  Voronoi  cells  associated  with  aggregate  elements  play  the  role  of  picture  elements, 
or  pixels.  Important  special  cases  include  hexagonal  cells  for  D  =  2,  and  truncated 
octahedral  cells  for  D  —  3.  The  cells  of  a  hexagonal  aggregate  are  centered  on  a 
hexagonal  lattice  and  the  cells  of  a  truncated  octahedral  aggregate  are  centered  on  a 
body-centered  cubic  lattice.  The  simplest  (level  1)  hexagonal  aggregate  can  be  visual- 
ized as  a  hexagonal  cell  with  six  neighboring  cells,  while  a  level  2  hexagonal  aggregate 
consists  of  49  cells.  A  central  cell  in  a  truncated  octahedral  aggregate  has  fourteen 
neighbors.  See  the  figures  in  Chapter  4.  The  cells  in  a  A-level  permutohedral  aggre- 
gate can  be  addressed  using  strings  of  length  A  from  the  set  of  digits  {0, 1, . . . ,  n  —  1}. 
This  form  of  addressing,  or  indexing,  is  known  as  generalized  balanced  ternary  (GBT), 
and  was  developed  by  Gibson  and  Lucas  [8]. 

The  majority  of  image  processing  techniques  and  algorithms  can  be  broadly 
classified  as  either  spatial  domain  or  frequency  domain  techniques.  Spatial  domain 
techniques  include  morphological  transforms,  while  frequency  domain  techniques  in- 
clude highpass  and  lowpass  filtering  (see  Gonzalez  and  Woods  [9],  or  Ritter  and 
Wilson  [18]).  Frequency  domain  techniques  depend  on  the  discrete  Fourier  trans- 
form (DFT),  a  linear  transform  that  extracts  the  frequency  content  of  images  with 
discrete  domains.  This  dissertation  is  concerned  with  the  computation  of  the  dis- 
crete Fourier  transform  and  the  inverse  DFT  of  functions,  or  images,  whose  domains 
are  permutohedral  aggregates.  The  geometric  duals  of  permutohedral  aggregates  are 
function  domains  in  the  frequency  space. 

The  fast  Fourier  transform  (FFT)  is  an  efficient  algorithm  for  the  computation 
of  the  discrete  Fourier  transform.  Gilbert  Strang  [19]  calls  the  FFT  "the  most  valuable 


numerical  algorithm  in  our  lifetime."  The  complexity  of  the  FFT  is  O(MlogM), 
where  M  is  the  cardinality  of  the  discrete  image  domain.  Mersereau  and  Speak  [15], 
and  Dudgeon  and  Mersereau  [7]  describe  the  FFT  for  periodic  functions  denned  on 
multidimensional  lattices,  and  Tolimieri,  An,  and  Lu  [20,  21],  describe  the  FFT  for 
functions  defined  over  Zn,  the  ring  of  integers  modulo  n.  The  FFT  is  based  on  the 
Cooley-Tukey  factorization  of  the  discrete  Fourier  transform,  which  was  published  in 
1965  [5].  However,  some  unpublished  manuscripts  by  Gauss  contained  the  essentials 
of  the  Cooley-Tukey  factorization  and  the  FFT  algorithm  around  1805  (see  Heideman 
et  al.  [11]). 

Data  reordering,  also  known  as  shuffle  permutations,  is  generally  associated 
with  FFT  algorithms.  Zhu  and  Ritter  [24,  26],  derived  algorithms  for  the  computa- 
tion of  the  FFT  of  images  defined  on  rectangular  arrays  that  did  not  require  a  data 
reordering  stage.  They  accomplished  this  by  using  the  generalized  matrix  product,  a 
generalization  of  standard  matrix  multiplication.  This  product,  also  known  as  the 
p-product,  was  first  described  by  Ritter  [17].  The  present  work  defines  the  gener- 
alized matrix  product  for  arrays  of  data  indexed  by  the  GBT  addressing  scheme. 
Such  a  definition  is  necessary  for  the  description  of  image  processing  algorithms  for 
permutohedral  aggregates  and  their  duals. 

This  dissertation  exploits  the  hierarchical  structure  of  permutohedral  aggre- 
gates and  the  GBT  addressing  scheme  to  derive  pyramidal  data  processing  algorithms 
that  compute  the  FFT  and  the  inverse  FFT  for  aggregates.  These  algorithms  have 
complexity  0(M  log  hi),  where  M  =  nx  and  n  =  2D+1  -  1;  they  are  expressed  in 
terms  of  the  generalized  matrix  product  for  permutohedral  aggregates. 


1.1     Hexagonal  and  Truncated  Octahedral  Aggregates 

We  include  some  further  remarks  regarding  hexagonal  and  truncated  octahe- 
dral aggregates  since  properties  of  these  two  special  cases  of  permutohedral  aggregates 
motivated  the  present  study. 

Watson  and  Ahumada,  Jr  [23]  derived  a  mathematical  transform  that  models 
space-frequency  coding  of  images  in  the  primary  visual  cortex  of  the  brain.  They 
noted  that  while  retinal  ganglion  cells  represent  visual  images  with  a  spatial  code,  the 
cells  of  the  primary  visual  cortex  code  spatial,  frequency,  and  orientation  information 
about  an  image.  Furthermore,  ganglion  cell  receptive  fields  form  an  approximate 
hexagonal  lattice  near  the  fovea,  the  region  in  the  central  portion  of  the  retina  that 
allows  human  beings  to  resolve  fine  details  in  an  image.  The  input  lattice  for  the 
image  transform  presented  by  Watson  and  Ahumada,  Jr.  is  a  hexagonal  array  that 
is  equivalent  to  the  hexagonal  aggregates  described  here. 

Petersen  and  Middleton  [16]  provide  further  motivation  for  the  study  of  multi- 
dimensional functions  sampled  on  hexagonal  and  body-centered  cubic  lattices.  They 
showed  that  the  most  efficient  sampling  lattices  are  not  in  general  rectangular.  Ef- 
ficient sampling  lattices  attempt  to  use  a  minimum  number  of  sampling  points  to 
exactly  reconstruct  bandlimited  functions.  In  particular,  sampling  two-dimensional 
isotropic  functions  on  a  hexagonal  lattice  and  sampling  three-dimensional  isotropic 
functions  on  a  body-centered  cubic  lattice  is  significantly  more  efficient  than  sam- 
pling such  functions  on  simple  square  and  cubic  lattices,  respectively.  The  defining 
property  of  an  isotropic  function  is  that  the  Fourier  transform  of  the  function  cuts  off 
at  the  same  magnitude  in  all  directions.  Thus,  the  support  of  the  Fourier  transform 
of  an  isotropic  function  is  a  hypersphere  or  ball.  Petersen  and  Middleton  specify 
minimum  sampling  lattices  for  isotropic  functions  in  dimensions  one  through  eight. 


1.2     Overview 

Chapter  2  presents  a  number  of  standard  definitions  and  results  from  finite 
abelian  group  harmonic  analysis:  the  discrete  Fourier  transform,  the  inverse  DFT, 
the  inversion  theorem,  the  convolution  theorem,  and  the  Cooley-Tukey  factorization 
theorem. 

Chapter  3  covers  the  discrete  Fourier  transform,  the  inverse  DFT,  and  the 
Cooley-Tukey  factorizations  for  finite  abelian  groups  that  arise  as  quotient  lattices 
and  their  geometric  duals,  respectively. 

Permutohedral  aggregates  and  the  quotient  lattices  A\  associated  with  per- 
mutohedral  aggregates  are  described  in  Chapter  4.  Permutohedral  aggregates  are 
particular  sets  of  coset  representatives  of  the  quotients  A\.  Chapter  4  also  discusses 
the  subgroups  Bp  in  A\,  which  lead  to  quotients  of  the  form  Qa  =  A\/Bp,  and  con- 
tains statements  of  the  Cooley-Tukey  factorizations  for  the  quotient  lattices  A\  and 
the  geometric  duals  A'x,  respectively. 

Chapter  5  defines  Y -indices — lexicographically  ordered  finite  strings  of  digits — 
and  the  generalized  matrix  product  over  T-indices.  The  T-indices  are  used  to  address 
the  elements  of  permutohedral  aggregates  and  their  duals. 

Two  theorems  in  Chapter  6  are  formulations  of  the  DFT  and  the  inverse  DFT 
for  permutohedral  aggregates  in  terms  of  the  generalized  matrix  product  over  T- 
indices,  or  digital  indexing  sets.  Chapter  6  also  states  four  main  digital  algorithms  in 
terms  of  the  generalized  matrix  product  for  permutohedral  aggregates:  the  Cooley- 
Tukey  (CT)  factorization,  the  inverse  CT  factorization,  the  fast  Fourier  transform, 
and  the  inverse  FFT. 

Chapter  7  contains  suggestions  for  further  research. 


1.3     Summary  of  Results 
The  main  results  established  in  this  work  are  as  follows: 

1.  algebraic  formulation  of  the  discrete  Fourier  transform  (DFT)  and  the  inverse 
DFT  for  permutohedral  aggregates; 

2.  definition  of  the  generalized  matrix  product  for  permutohedral  aggregates; 

3.  discrete  Fourier  transform  and  inverse  DFT  algorithms  in  terms  of  the  gener- 
alized matrix  product  for  permutohedral  aggregates; 

4.  Cooley-Tukey  factorization  and  inverse  CT  factorization  algorithms  in  terms  of 
the  generalized  matrix  product  for  permutohedral  aggregates; 

5.  a  decimation-in-space,  fast  Fourier  transform  (FFT)  algorithm  for  permutohe- 
dral aggregates; 

6.  a  decimation-in-frequency,  inverse  FFT  algorithm  for  permutohedral  aggre- 
gates. 


CHAPTER  2 
THE  DISCRETE  FOURIER  TRANSFORM 


2.1     Introduction 

The  purpose  of  this  chapter  is  to  present  background  material  from  finite 
abelian  group  harmonic  analysis.  Section  2.2  describes  the  normed  vector  space  L2(A) 
of  complex-valued  functions  on  a  finite  abelian  group  A,  the  algebraic  dual  A,  and 
the  space  L2(A).  Section  2.3  presents  some  standard  definitions  and  results  for  the 
discrete  Fourier  transform  and  the  inverse  DFT  for  finite  abelian  groups.  We  recall 
that  the  Fourier  transform  is  an  isometry  and  state  the  inversion  theorem.  Section  2.4 
defines  cyclic  convolution  for  functions  in  L(A)  and  presents  the  convolution  theorem 
for  finite  abelian  groups.  In  Section  2.5  we  prove  the  Cooley-Tukey  factorization 
theorem. 

Denote  the  set  of  positive  integers  by  Z+  and  the  field  of  complex  numbers 
byC. 

2.2     Vector  Spaces  Over  Finite  Abelian  Groups 

Denote  by  card(5)  the  cardinality  of  a  set  S.  Let  A  be  a  finite  abelian  group. 
The  set  L(A)  of  all  complex-valued  functions  defined  on  A  is  a  vector  space  over  C 
with  addition  and  multiplication  given  by 

(a  +  b)(r)  =  a(r)  +  6(r),  for  all  a,  6  €  L(A),  r  e  A, 
(ta)(r)  =  £a{r),  for  all  a  6  L(A),  £  6  C,r  €  A. 


For  a, be  L(A),  define  ab  €  L(A)  by  (ab)(r)  =  a(r)b(r),  for  all  r  €  A. 
Define  an  inner  product  on  L(A)  as  follows.  For  a,  b  €  £(-4),  let 

{a,6)  =  £a(r)6B 


where  6(r)  denotes  the  complex  conjugate  of  6(r).  For  a  €  £(>!),  define  the  norm  of 
a  by 


\\a\\  =  y^)  =  (Y\a(r)A 


1/2 


Denote  by  £2(j4)  the  inner  product  space  L(A)  together  with  the  norm  ||  •  |j. 

We  now  describe  a  basis  of  L2(A).   For  r,  s  £  A,  define  the  Kronecker  delta 

over  indices  in  A  by 

r         J  1    if  r  =  s 

r°  ~  \  0    if  r  ^  s. 

For  r  £  A,  define  er  6  //(>!)  by  er(s)  =  5r3,  for  all  s  E  A.  We  will  also  have  occasion 
to  refer  to  the  support  of  a  function,  defined  as  follows.  If  /  is  any  function  with 
domain  A,  R  is  any  ring,  and  /  :  A  — >  R,  then  the  support  of  /,  written  supp(/),  is 
defined  by  supp(/)  =  {d  6  A  :  f(d)  ^  0}.  Thus,  supp(er)  =  r,  for  all  r  6  A. 

Proposition  1   The  set  S  =  {er  :  r  €  A]  is  an  orthonormal  basis  of  L2(A). 
Proof.     Ha£L2(A),  then 


a  =  \Ja(r)er. 


re.4 

Hence,  the  set  of  vectors  {er  :  r  €  ^4}  spans  the  space  L2(A). 
If  r  and  s  are  distinct  elements  of  A,  then 

(er,  e4)  =  ^  er(t)es{t)  =  er(s)  +  e,{r)  =  £.,  +  £4r  =  0. 

Thus,  £  is  a  set  of  (nonzero)  orthogonal  vectors,  which  implies  that  it  is  a  linearly 
independent  set.  Since  £  also  spans  L2(A),  it  is  a  basis  of  L2(A). 


If  r  —  s,  then 

(Cp,e,)  =  (er,er)  =  ^2  er(t)er(t)  =  1. 

Hence,  £  is  an  orthonormal  basis  of  L2(A).  D 


Refer  to  £  as  the  standard  basis  of  £2(A).  If  card(A)  =  n,  then  £2(j4)  is  an 
n-dimensional  normed  vector  space,  since  card(£)  =  card(yl).  Section  2.3  describes 
a  second  basis  of  L2(A). 

Let  U  =  {z  6  C  :  \z\  =  1},  the  unit  circle,  and  for  n  €  Z+,  let  f/(n)  denote 
the  multiplicative  group  of  all  the  nth  roots  of  unity.  Define  the  algebraic  dual  A 
of  a  finite  abelian  group  A  of  cardinally  n  by  A  =  Hom(  A  U(n)),  the  set  of  all 
homomorphisms  from  A  into  U(n).  Thus,  (f)  £  A  implies  (f>  :  A  — >  U(n)  such  that 
(f>(r  +  s)  =  <f>(r)(f)(s),  for  all  r,s  £  A.  Define  addition  in  A  as  follows:  for  d>,  rp  €  A,  let 
(0  +  0)(r)  =  <£(r)</>(r),  for  all  r  6  A.  Note  that  k<t>{r)  =  4>(r)k,  for  all  A:  €  Z,r  €  /I. 
This  addition  renders  ^4  an  abelian  group.  Denote  by  0o  the  zero  element  or  trivial 
homomorphism  in  A  Hence,  <£o(r)  =  1,  for  all  r  £  A. 

A  standard  result  for  finite  abelian  groups  follows. 

Theorem  2  If  A  is  a  finite  abelian  group  with  card(A)  =  n  and  A  =  Hom(A,  U(n)), 
then  A  =  A.  In  particular,  card(A)  =  n. 

The  elements  of  the  dual  group  A  enjoy  the  following  properties. 

Proposition  3  If  A  is  a  finite  abelian  group  and  4>  €  A,  then 

1.  0(0)  =  1; 


2.   \<f>{r)\2  =  <f>{r)4>{r)  =  1,  for  all  r  €  A; 


3.   The  additive  inverse  —<j)  of  <f>  satisfies  —<p(r)  =  <f>(r)   1  =  4>(r),  for  all  r  €  A 
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Proof.  Any  homomorphism  <p  :  A  ->  U(n)  maps  the  identity  element  0  €  A  to  1, 
the  identity  element  of  £/(n).  Hence,  <£(())  =  l,for  all  <f>  €  A.  which  proves  (1). 

If  r  e  A,  then  |0(r)|2  is  the  norm  of  <j>(r)  €  U(n).  The  elements  of  U{n)  lie 
on  the  unit  circle.  Hence,  |</>(r)|2  =  l,for  all  r  €  A,  which  proves  (2). 

If  —  <p  is  the  additive  inverse  of  <f>  €  A,  then  (<£  +  (— <£))(r)  =  <£o(r)  = 
l,for  all  r  6  A.    Since  (dt  +  (-<£))(r)  =  (0(r))(-0(r))  =  1,  -0(r)  is  the  multi- 


plicative inverse  of  <p(r)  €  £/(n).   Hence,  —  0(r)  =  (/>(r)   '  and  <f)(r)   x  =  <£(r),  for  all 
r  €  A,  since  —  <£(r)  €  £/,  for  all  r  €  A,  which  proves  (3).  D 

Note  that  <j>  €  A  implies  <f> :  A  -*  U(n)  C  C,  so  that  d>  €  L(A),  for  all  <£  €  A 
Thus,  Ac  L(A). 

The  set  L(A)  =  {a  :  A  — >  C}  of  all  complex- valued  functions  on  /I  is  a  vector 
space  over  C  with  addition  and  scalar  multiplication  given  by 

(a  +  b)(d>)  m  a(d>)  +  b{<f>),  for  all  a,  6  €  £(A),  </>  €  A, 
((S)(<f>)  =  (8{<t>),  for  all  a  €  1(A),  £  €  C,  <j>  €  A. 

For  a,  b  €  L(A),  define  a  •  6  €  L(A)  by  (a  ■  6)(0)  =  a{4>)b{<p),  for  all  <£  €  A. 
Define  an  inner  product  and  a  norm  on  L(A)  as  follows.  For  o,i  £  A,  the  inner 
product  of  a  and  6  is  given  by 


(2,6)  =  5>(#(<£), 
and  the  norm  of  a  is  given  by 

iio||  =  v^)=  (2>(*)ia)    • 

In  Section  2.3  we  define  a  basis  of  L(A)  with  cardinality  n.  Denote  by  L2(A) 
the  inner  product  space  L(A)  together  with  the  norm  ||  •  ||.  Thus,  L?{A)  is  a  normed 
vector  space  over  C  of  dimension  n. 
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In  order  to  define  the  Fourier  transform  for  finite  abelian  groups  in  Section  2.3 
we  need  a  bilinear  form  (•,•):  A  x  A  -t  U(n). 

Proposition  4  If  A  is  a  finite  abelian  group  of  cardinality  n,  then  the  map 

{•}-):AxA—+U(n) 

defined  by 

(r,4>)  =  <}>{r)i 

for  all  r  €  A  and  d>  6  A,  is  a  bilinear  form,  that  is, 

1.  (r  +  5,0)  =  (r, <f>)(s, 4>),  for  all  r,s  €  A,d>£  A; 

2.  (r,  4>  +  V>)  =  (r,  <j>)(s,  <j>),  for  all  r  e  A,<t>,xp  6  A; 

3.  (kr,  <j>)  =  (r,  <f>)k  =  (r,  kd>),  for  all  r  €  A,d>€  A. 


In  particular,  Proposition  4  implies  (—r,<f>)  =   (r, d>)       =  (r, <f>)  =  (r,—<f>), 
for  all  r  €  A,  d>  €  A. 

2.3     The  Discrete  Fourier  Transform  and  Inverse  DFT 

Definition  1  (Discrete  Fourier  Transform)  If  A  is  a  finite  abelian  group  with 
card(A)  =  n,  then  the  Fourier  transform  of  a  €  L(A)  is  the  function  a  6  L(A) 
defined  by 

S(*)  =  4jX>(0(r,*),  (2-1) 

/or  all  (j>  6  A.   W^e  s/ia//  tvn'te  ^  /or  f/ie  map 

^  :  1(A)  — *  L(A) 
such  that  Ta{o.)  =  a. 
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The  complex  numbers  a{4>)  are  known  as  the  Fourier  coefficients  of  a  6  L(A). 
Computation  of  the  DFT  of  a  €  L(A)  via  equation  2.1  requires  no  more  than  n2 
complex  multiplications  (not  counting  multiplication  by  the  normalization  constant). 
The  factors  (r,  0)  can  be  kept  in  a  lookup  table.  Hence,  their  computation  does  not 
contribute  to  the  complexity  of  equation  2.1. 

Example.  For  a  €  L(A)  and  s  €  A,  the  shift  of  a  by  s  is  the  function  a,  6  L(A) 
defined  by  a„(r)  =  a(r  —  s),  for  all  r  €  A.  Applying  the  Fourier  transform  to  a,  yields 
as{4>)  =  {s,<fr)  -a{(t>)  =  <p(s)  ■  a(<f)),  for  all  0  €  A. 

One  may  easily  verify  linearity  of  the  Fourier  transform,  which  we  state  for- 
mally. 

Proposition  5  If  A  is  a  finite  abelian  group  with  card{A)  =  n,  then  Ta  is  a  linear 
transform,  that  is, 

7A{a  +  b)  =  TA(a)  +  FA{b),  for  all  a,be  L{A), 

and 

FA(t*)  =  6^(a),  for  all  a  €  L(A),£  €  C. 

For  an  application  of  the  Fourier  transform,  we  determine  a  basis  of  L2(A) 
by  computing  the  Fourier  transforms  of  the  basis  elements  in  £  =  {er  :  r  €  A}.  Let 
er  =  FA(er),  for  all  r  €  A. 


Lemma  6  If  r  £  A,  then 


for  all  <j>  €  A. 


eM)  =  -7=4>{r), 
\/n 
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Proof.    If  r  €  A  then 


.{cj>)    =    -L£er(5)(5,0) 

V  mC   A 


for  all  0  €  A.  □ 


Noteworthy  is  the  function  e0,  where  0  is  the  zero  element  of  A.  Since 

e0(4>)  =  ±m  =  -L, 

y/n  y/n 

for  all  (f)  G  A,  e"o  is  constant  on  A. 

Lemma  7  If  r  £  A,  then 

E(r  ±\  _  J  n-    i/>  =  0 
V  ,W       10,     otherwise. 


4>eA 
Proof.     If  r  =  0,  then  (r,  </>)  =  (0,<£)  =  ^(0)  =  1,  for  all  <j>  €  A   The  result  follows 
from  the  cardinality  of  A.  If  r  ^  0,  then  there  exists  tp  £  A  such  that  ip(r)  ^  1,  that 
is,  (r,  V0  ^  0.  Thus, 


4>€-4 
=      (r»^)£(r>^)l 


^ev4  </>€,4 


*€>! 


which  implies 


<tf>€-4 
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Proposition  8    The  set  £  =  {er  :  r  £  A}  is  an  orthonormal  basis  of  L2(A). 

Proof.     Suppose  r,s  €.  A.   Using  the  inner  product  in  L2(A)  and  the  bilinear  form 
on  A  x  A,  we  have 

(er,es)    =    ^er(<£)e,(0) 
^e.4 

=     —  y    <f>(r)4>(s)         (by  Lemma  6) 
n  ^— ' 

4>£A 

06^4 

77   * — ' 

1,    if  r  =  5 
0,    otherwise, 

the  last  equality  by  Lemma  7.  Thus,  (er,e3)  =  <£r4,  for  all  r,  5  6  A,  which  implies  £" 

is  an  orthonormal  basis  of  L2(A).  D 

Since  the  cardinality  of  £  is  n,  L(A)  has  dimension  n.  Furthermore,  since  Ta 
maps  an  orthonormal  basis  to  an  orthonormal  basis,  we  have  the  following  result. 

Theorem  9   The  Fourier  transform  Ta  '•  L2(A)  — ►  L2(A)  is  an  isometry,  that  is, 

\\FA{a)\\  =  \\a\\,foralla£L2{A). 

Before  defining  the  inverse  Fourier  transform,  we  present  a  dual  to  Lemma  7. 
Lemma  10  If  <f>  €  A,  then 


rg.4 

where  <f>0  is  the  trivial  homomorphism 


\f  <j>  =  (f>Q 

therwise, 
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Proof.    If  <j>  =  <£o,  then  (r,  <£)  =  4>o{r)  =  1,  for  all  r  €  A  The  result  follows  from  the 
cardinality  of  A.  If  0  ^  </>0,  then  <£(s)  ^  1,  for  some  s  6  A,  and 

reA 


r€>4 


implies 


Dr^)  =  °- 


r€>4 


The  next  example  yields  another  basis  of  L2(A). 
Example.  If  0,  ip  €  A,  then  4>,ip  £  L(A).  We  compute  the  inner  product  of  (f>  and 


(*,V)     =     £>(r)V(r) 

r£A 


r£A 

n,    if  (/>  =  ip 
0,     otherwise, 

the  last  equality  by  Lemma  10.  Normalization  yields  an  orthonormal  basis 


€'  =  {-L/> :  <f>  e  A}  =  -±=A 


of  L2{A). 
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Example.  If  4>  €  A,  then  4>  £  L(A).  What  are  the  Fourier  coefficients  of  <f>  =  FaW*- 
If  tl>  €  A,  then 

J(0)    =     -i=5>(r)(r,0) 
V  n  'T^ 

r£A 

=    -f=y2(r,<J>)(r,ti>) 

f  Vn.    if  V  =  — ^ 
|  0,        otherwise, 

the  last  equality  by  Lemma  10. 


Definition  2  (Inverse  Discrete  Fourier  Transform)  If  A  is  a  finite  abelian  group 
with  card(A)  =  n,  A  is  the  algebraic  dual  of  A,  and  a  6  L(A),  then  the  inverse  Fourier 
transform  of  a  is  the  function  a  =  J~2  (<*)  €  L{A)  given  by 


yjn  *— i 

4>£A 


for  all  r  E  A. 


Computation  of  the  inverse  DFT  of  a  €  L(A)  via  equation  2.2  requires  at 
most  n2  complex  multiplications. 

The  composition  of  Ta  and  T~^  is  the  identity  transformation  on  L(A),  as 
stated  in  the  next  result. 

Theorem  11  (Inversion  Theorem)  If  a  €  L(A)  and  a  =  Ta{<"),  then  a  =  ^'(a). 
Proof.     Let  a  €  L(A)  and  b  =  ^(a).  If  s  €  A,  then 

06-4 
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=    iE(£S>Xr,#))E« 

4>eA    \v        rCA  / 

n  * — '  * — ' 

4>€Ar^A 

=     —a{s)  •  n         (by  Lemma  7) 
n 

=    a(s), 
which  implies  b  =  a.  □ 


2.4     The  Convolution  Theorem 

Definition  3  (Convolution)  If  a,b  €  £(>!),  #ien  the  (cyclic)  convolution  c  €  L(A) 
of  a  and  b,  written  c  =  a  *  b,  is  defined  by 


c(0  =  £a(r)6(<-r), 

reA 


for  all  t  e  A. 


Example.   We  compute  the  convolution  of  a  basis  element  es  €  €  with  a  function 
a  €  L(A)  at  t  €  A: 

(e,  *  a)(t)  =  E  e*(r)a(t  -  r)  =  a(t  -  s)  =  as(t), 

reA 

where  a3  is  the  shift  of  a  by  s  €  A.  Thus,  es  *  a  =  at. 
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Example.  If  <j>,  tp  €  A,  then  0,0  €  £(4)  and 

(0*0)(O)    =    £>(rMO-r) 

=      $^(r><£)(-ri0) 


r6^ 

n,    if  ^  =  %l> 
0,     otherwise, 


the  last  equality  by  Lemma  10. 


The  next  example  shows  that  we  may  regard  the  Fourier  transform  as  a  con- 
volution. 

Example.  If  a  €  L(A)  and  <f>  €  A,  then  the  properties  of  the  bilinear  form  on  A  x  A 
imply 


Thus, 


(a*(-0))(O)    =    5>(r)(-#0  -  r)) 

reA 

=    ^2a(r){-r,-<l>) 

reA 

=    $3«(r)(r.^) 

reA 

=     y/n  ■  a((j>). 


a(<f>)  =  ^=-{a*  (-*))(0). 


The  Fourier  transform  relates  convolution  and  pointwise  multiplication.  This 
relation  is  known  as  the  convolution  theorem. 


19 

Theorem  12  (Convolution  Theorem)  If  A  is  an  abelian  group  with  card(A)  =  n, 
a,  6  €  L2{A)  and  c  =  a*  b,  then  c  =  y/n  -ah. 

Proof.     If  a,b,c  G  L(A),  with  c  =  a  *b,  and  <f>  €  A,  then 

(a-b)(4>)     =    a(4>).b(4>) 

^a(r)(r,^^^fe(5)(,» 

reA  \seA  / 

r€A    \s£A  J 

n  L — '  L — ' 
reA  seA 

=  -EEa(r)6(s)(5+r^)' 

n  * — »  ' — • 
using  bilinearity.  Letting  £  =  s  +  r  and  interchanging  the  summations, 

(a.b)(4>)    =    i]Ty>(r)6(*-r)(*,*) 

reA  teA 

=  i  £<*<)(*,*) 

r)   '     • 


n 

teA 

Vn 


Example.  Since  as  =  ea  *  a  from  a  previous  example,  the  Fourier  transform  of  a„  is 
a*  =  y/n  •  e,  •  a,  by  the  convolution  theorem. 
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2.5     The  Cooley-Tukey  Factorization 

The  Cooley-Tukey  (CT)  factorization  is  an  algorithm  that  computes  the  dis- 
crete Fourier  transform  of  a  function  a  6  L(A),  where  A  is  a  finite  abelian  group.  An 
immediate  corollary  computes  the  inverse  DFT  of  a  function  a  €  L(A). 

For  the  spatial  domain,  consider  the  short  exact  sequence 

0  — ►  B  — ►  A  — ->  A/B  — ►  0, 

where  A  is  a  finite  abelian  group  and  B  is  a  subgroup  of  A.  Note  that  B  does  not 
necessarily  split  in  A,  that  is,  it  is  not  generally  true  that  there  exists  a  subgroup  C 
of  A  such  that  A  =  B  ©  C  (internal  direct  sum).  However,  if  R  C  A  is  a  complete 
and  nonredundant  set  of  coset  representatives  of  the  quotient  group  A/B  such  that 
B  D  R  =  {0},  then  A  =  B  +  R,  and  for  every  r  €  A  there  exist  unique  rp  €  B  and 
rQ  €  R  such  that  r  —  rp  +  ra.  We  refer  to  such  a  set  as  a  complete  system  of  coset 
representatives  (or  residues)  of  A/B. 

Let  Bm  =  {4>  £  A  :  (r,  0)  =  1,  for  all  r  €  B}.  It  is  easy  to  see  that  B»  is  a 
subgroup  of  A.  In  the  frequency  domain  consider  the  short  exact  sequence 


4>o  — ►  5,  — ►  A  — ►  A/5,  — ►  $ 


'o> 


where  <£0  is  the  trivial  homomorphism  in  A.  If  i2,  C  A  is  a  complete  system  of 
coset  representatives  of  A/B,,  then  A  =  Bm  +  R,  and  for  every  <f>  G  A  there  exist 
unique  0Q  e  5,  and  ^  €  R,  such  that  0  =  (j)a  +  d>0.  One  can  show  that  A/Bm  =  B 
and  £.  £  A/5.  Let  n  =  card(A),  n/j  =  card(£),  and  na  =  card(A/5).  We  have 
n  =  nan0  =  card(A),  card(R)  =  na  =  card(5»),  and  card(A/5.)  =  nfi  =  card(#„). 
An  et  al.  [1],  give  a  detailed  discussion  of  the  next  result,  although  their  approach 
depends  on  an  isomorphism  from  A  onto  A. 
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Theorem  13  (Cooley-Tukey  Factorization)  Let  A  be  a  finite  abelian  group,  B 
a  subgroup  of  A,  and  R  C  A  a  complete  system  of  coset  representatives  of  A/B.  If 
a  €  L(A)  and  a  =  J-a{u),  then 

S(*)  =  4=£   f  f-^£a(r)(^,^)]-(r0,^)]  -(ra,0a),  (2.3) 

for  all  4>  =  d>a  +  dtp  €  B»  +  R*  =  A,  where  R*  C  A  is  a  complete  system  of  coset 
representatives  of  A/B*,  and  r  =  r$  +  ra  £  B  +  R  =  A. 

Proof.  Let  d>  =  4>a  +  <pp  €  A  with  4>a  G  B.  and  <£/?  €  /?..  If  r#  G  B,  then  (r^,  <fc)  =  1 
by  the  definition  of  B..  Hence, 

(r,4>)    =    {r0  +  ra,dia  +  <f>p) 

=    {r0,d>a)(r0,d>0)(ra,4>a)(ra,<f>p) 
=    (r0,<t>0)(ra,d>0)(ra,d>a), 

for  all  r  =  rp  +  rQ  €  A  and  <f>  =  4>a  +  d>0  €  A.  The  Fourier  transform 

3W--jffE-MM).  (2-4) 

where  n  =  n0na  is  the  cardinality  of  A,  can  therefore  be  written  as  in  the  statement 
of  the  theorem  (equation  2.3).  □ 


Corollary  14  Let  A  be  a  finite  abelian  group,  B  a  subgroup  of  A,  and  R,  C  A  a 
complete  system  of  coset  representatives  of  A/B*.  If  a  €  L(A)  and  a  —  ^^(a),  then 

a{r)  =  —    ]T    I  I  —    Y^   a{4>){r*,4>a)  )  ■  (ra,+fi)  I  ■  (r^,^), 

/or  all  r  =  r0  +  ra  €  B  +  i?  =  .4,  w/iere  RC  Au<  complete  system  of  coset 
representatives  of  A/B,  and  (f)  =  4>a  +  <f>0  €.  B*  +  R,  =  A. 
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The  numbers  (ra,^)  and  {ra,<f>p)  in  Theorem  13  and  Corollary  14,  respec- 
tively, are  known  as  twiddle  factors.  The  Cooley-Tukey  factorization  is  a  three-stage 
algorithm.  The  first  stage  is  the  inner  summation  in  equation  2.3  which  requires  no 
more  than  npna  =  n  complex  multiplications  for  each  <j>p,  for  a  total  of  n  ■  np  complex 
multiplications.  The  second  stage  is  an  entry-by-entry  multiplication  by  the  twiddle 
factors  consisting  of  nanp  complex  multiplications.  The  last  stage  is  the  outer  sum- 
mation in  equation  2.3  which  requires  no  more  than  nanp  complex  multiplications  for 
each  <j>a,  for  a  total  of  n  ■  na  complex  multiplications.  The  total  number  of  complex 
multiplications  is  bounded  by  n  ■  np  +  n  ■  na  +  nanp.  Consider  the  computational 
savings  for  the  case  n  —  104,np  =  100  =  na.  Direct  computation  using  equation  2.4 
involves  n2  =  108  complex  multiplications.  The  Cooley-Tukey  factorization  algorithm 
requires  no  more  that  2  •  106  +  104  complex  multiplications. 


CHAPTER  3 
QUOTIENT  LATTICES 


3.1     Introduction 


In  this  chapter  we  consider  the  discrete  Fourier  transform  for  functions  defined 
on  finite  abelian  groups  that  arise  as  quotients  of  /^-dimensional  lattices.  Instead  of 
working  with  algebraic  duals  in  the  frequency  domain,  we  will  take  a  more  intuitive 
approach  and  work  with  the  equivalent  geometric  duals.  This  approach  differs  from 
standard  presentations  of  the  discrete  Fourier  transform  that  emphasize  the  algebraic 
structure  of  the  underlying  indexing  sets  for  function  domains;  our  emphasis  is  on 
the  geometric  structure  of  the  underlying  indexing  sets  or  coset  representatives. 

Section  3.2  presents  the  discrete  Fourier  transform  and  the  inverse  DFT  for 
quotient  lattices.  Since  quotient  lattices  have  the  structure  of  finite  abelian  groups, 
all  the  results  of  Chapter  2  are  immediately  applicable.  In  particular,  the  inversion 
and  convolution  theorems  apply  to  the  quotients  under  consideration  in  this  chapter. 
Section  3.3  states  the  Cooley-Tukey  factorization  for  quotient  lattices  in  the  spatial 
and  frequency  domains. 

Denote  the  field  of  real  numbers  by  R,  Z)-dimensional  Euclidean  space  by  RD, 
and  the  set  of  integers  by  Z.  Denote  the  inner  product  of  x,  y  E  RD  by  (x,y). 

3.2     The  DFT  and  Inverse  DFT  for  Quotient  Lattices 

A  /^-dimensional  lattice  is  the  set  of  all  linear  integer  combinations  of  D 
linearly  independent  vectors  in  R    .  Let  L  be  a  Z)-dimensional  lattice.  The  geometric 
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dual  (or  reciprocal  lattice)  L*  of  L  is  defined  by 

L'  =  {s  <E  RD  :  (r,  3)  €  Z,  for  all  r  €  L}. 

Note  that  Z*  is  also  a  /^-dimensional  lattice.  Denote  a  sublattice  relationship  by 
L\  <  Lq  for  two  Z)-dimensional  lattices  L0  and  L\.  Note  that  "-<"  is  also  an  abelian 
subgroup  relationship. 

Proposition  15  (Lattice  Duality  Relationship)  If  Lq  and  Li  are  D -dimensional 
latttices,  then 

L\  -<  Lq  «(     V  Lq  -K  Lj 

Remark  1  From  a  signal  processing  perspective,  we  regard  Lq  in  Proposition  15  as 
a  spatial  domain  sampling  lattice  and  L\  as  a  spatial  domain  periodicity  lattice. 
Similarly,  we  regard  L\  as  a  frequency  domain  sampling  lattice  and  Lq  as  a  frequency 
domain  periodicity  lattice. 

The  sublattice  relationship  also  applies  to  sublattices  of  quotient  lattices  of 
the  form  LqJL\.  Quotient  sublattices  arise  from  sequences  of  the  form  L\  <  L2  ■*<  L0. 

Corollary  16  If  Lq,  L\,  and  L2  are  D-dimensional  latttices  such  that  L\  -<  L2  -<  Lq, 

then 

1.      Lq    -~s    L/2    "*\    **j/ 

2.  L2/L1  -<  Lo/Li  and  L^/L*0  <  L*/L*0; 

3.  {LofLx)l{LifLx)  3  Lq/L2  and  (L*J L'q)/(L*2/L*q)  3  L'JL;. 

Define  quotient  lattices  for  the  spatial  and  frequency  domains  by  A  =  Lq/ L\ 
and  A'  =  L\/ Lq,  respectively.  If  n  is  the  index  of  Lt  in  L0,  then  A  is  a  finite  abelian 


25 

group  with  cardinality  n,  and  nr  €  Li,for  all  r  €  L0.  By  Proposition  18  (below),  the 
index  of  L*Q  in  L\  is  also  n.  Hence,  ns  6  Lo,for  all  5  €  LJ.  The  spatial  and  frequency 
domain  short  exact  sequences  corresponding  to  A  and  A'  are 

0  — >  Li  — >•  L0  — ►  A  — >0 

and 

0  — ►  Lo  — ►  L\  — >  A'  — ►  0, 

respectively. 

Recall  that  (7(n)  denotes  the  multiplicative  group  of  all  the  nth  roots  of  unity. 
The  next  result  defines  a  bilinear  form  on  A  x  A'. 

Proposition  17  Let  LQ  and  Lx  be  D -dimensional  lattices  with  L\  -<  L0.    If  A  = 

Lo/Li,  A'  =  Li/Lq,  and  n  =  card(A),  then  the  map 

{■,-):AxA'^>U{n) 


defined  by 

(r,s)  =  exp[— 27ri(r,  s)], 

where  i  =  \/—l  and  f  and  s  represent  the  cosets  containing  r  and  s,  respectively,  is 
a  bilinear  form. 

Proof.  Since  A  =  L0/Li  and  nr  e  L\,  for  all  r  6  Z-o,  we  have  (nr, s)  6  Z,  for  all 
s  E  L\,  that  is,  (r,s)  G  *Z,  for  all  r  6  L0  and  s  e  L\,  which  implies  (r,s)  = 
(r  +  Li,5  +  Lq)  =  exp[— 27ri(r,  5)]  G  £/(n).  To  show  that  (r,s)  is  independent  of  the 
choice  of  coset  representatives,  let  ri,r2  €  Lq  with  r\  =  r2  mod  L\.  Then  n  —  r2  =  u, 
for  some  u  £  L\.  Now,  if  5  G  ZJ,  then 

exp[-2ni(ri, s)}  =  exp[-27ri(r2,  s)]     «=^    exp[-27ri(r1  -  r2,  s)]  =  1 

<*=>•    exp[— 2iri(u,s)]  =  1, 
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which  is  true  for  all  s  €  L\  since  u  €  L\.    A  similar  argument  shows  that  (r,s)  is 
independent  of  the  choice  of  coset  representative  of  the  equivalence  class  s.    Thus, 
(•,  •)  is  well  defined. 
The  identities 

(n  +  r2,s)  =  (ri,s)(r2,s),      for  all  rur2  €  A,  s  €  A', 

(F,5!  +  s2)  =  (r,sx)(r,s2),      for  all  r  €  A,  Si,s2  €  A', 

and 

(cr,  s)  =  (r,  s)c  =  (r,  cs),      for  all  c  g  Z,  r  €  A,  s  €  A' 

are  all  easily  verified.  □ 

Note  that  (r,  5)  =  1  if  and  only  if  (r,  s)  €  Z,  where  r  is  any  coset  representative 
oir^A  and  s  is  any  coset  representative  of  s  €  A'. 

The  quotient  lattice  A'  is  the  geometric  analogue  of  the  dual  group  A  = 
Hom(/l,  U(n)).  Two  abelian  group  isomorphisms  follow. 

Proposition  18  7/  L0  and  L\  are  D -dimensional  lattices  such  that  L\  ■<  Lq,  A  = 
L0/Li,  A'  =  L\/Lq,  and  n  =  card(A)  (as  in  Proposition  17),  then  (1)  A'  =  A,  and 
(2)  A'  ^  A. 

Proof.  Part  (1).  For  every  s  -  s  +  L"Q  G  A',  define  a  map  fa  :  A  ->  U(n),  where 
n  =  card(A),  by  <^(r)  =  (r,s),  for  all  r  =  r  +  L\  €  A.  The  linearity  of  (-,s)  implies 
fa  £  A,  that  is,  fa  is  a  homomorphism.  Define  $  :  A'  —>  A  by  ty(s)  =  fa,  for  all 
5  €  A'.  The  linearity  of  (r,  •)  implies  ^  is  a  homomorphism.  Since  card(A)  =  n  = 
card(v4)  (Theorem  2),  we  need  only  show  ty  is  injective.  If  Sx,s2  €  A',  then 

tffa)  =  tf(s2)     =*     (r,li)  =  (r,s2),  for  all  r  €  A 
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=*  (r,3i)(r,32)_1  =  1 

=►  (r,6i)(r, -52)  =  1 

=*  (r,3i-32)  =  l 

==>•  exp[—2ni(r,si  —  s2)]  =  1,  for  all  r  €  L0, 

=>  Si  =  s2 

=►  3l   =  52. 

Hence,  $  is  an  isomorphism. 

Part  (2).  By  Theorem  2,  A  =  A;  by  the  first  part  of  the  proposition,  A'  =  A. 
Hence,  A'&  A.  D 

Denote  by  L(A)  the  n-dimensional  space  of  all  complex-valued  functions  on 
A  and  by  L(A')  the  n-dimensional  space  of  all  complex-valued  functions  on  A'.  The 
geometric  analogue  of  the  Fourier  transform  (Definition  1)  of  a  function  a  €  L(A)  is 
the  function  a'  =  Fa{o),  where  Ta  '■  L(A)  — >  L (A'),  defined  by 

a'(s)  =  -=  ^2a(r){r,s)  =  —=  ^a(r)exp[-27ri(r,«)],     for  all  s  €  A'.      (3.1) 

The  inverse  Fourier  transform  of  a  function  a'  €  L(A')  (the  geometric  analogue  of 
Definition  2)  is  the  function  a  =  ^jt  (<*')  €  £(AJ  defined  by 

o(r)  =  -7=  T]  o'(5)(r,I)  =  -7=  V  a'(s)  exp[2™'(r, «)],     for  all  r  €  A       (3.2) 

In  the  sequel  we  will  work  with  sets  of  coset  representatives  of  A  whenever  A 
is  a  quotient  lattice,  that  is,  whenever  A  =  L0/Li  for  Z)-dimensional  lattices  L0  and 
L\.  Similarly  for  the  dual  quotient  A'  =  L\/L$.  In  this  context,  let  G  C  L0  be  a 
complete  system  of  residues  for  A  and  let  G'  C  L\  be  a  complete  system  of  residues 
for  A'.  Assume  0  €  G  and  0  €  G",  where  0  is  the  zero  vector  in  RD.  Define  a  map 
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(•,•)  on  L0  x  L\  by 

(r, 5)  =  exp[-27n'(r, *)],     for  all  r  €  L0,s  €  LJ.  (3.3) 

Since  L0  and  Lj  are  Z-modules,  and  from  the  proof  of  Proposition  17,  we  have 

(v):Loxi;— ►*/(„), 

where  n  =  card(A),  is  a  bilinear  form.  Now,  define  a  :  G*  — ►  C  by  a(r)  =  a(r),  for  all 
r  <E  G\  and  a' :  G"  ->•  C  by  a'(s)  =  a'(s),  for  all  5  €  G".  Note  that  G  x  G'  C  L0  x  L\. 
In  view  of  Proposition  17,  compute  the  DFT  of  a  €  L(A)  as 

a'(s)  =  -i=^2  a(r)(r' s)  =  "7"  5Z  a(r^  exP[-27r*'<r>  *>li     for  a11  5  €  G",        (3.4) 
and  the  inverse  DFT  of  a'  6  £(,4')  as 

a(r)  =  —/=^2  a'(5)(r'5)  =  -7=  5Z  a'(5)  exP[27"'(r>  s>]>     for  a11  r  eG.        (3.5) 

3.3     The  Cooley-Tukey  Factorization  for  Quotient  Lattices 

The  purpose  of  this  section  is  to  state  the  Cooley-Tukey  (CT)  factorization 
(Theorem  13)  and  the  inverse  CT  factorization  (Corollary  14)  for  the  pair  of  sequences 
{B  — ►  A  — >■  A/B,  B'  — ►  A'  — ►  A'/B'},  respectively,  where  A  and  A1  arise  as  quotients 
of  D-dimensional  lattices,  B  is  a  subgroup  of  A,  and  B'  is  a  subgroup  of  A'.  If 
np  =  card(B)  and  na  =  ca,vd(A/B),  then  n  =  npria  =  card(A). 

In  the  spatial  domain,  let  Lq,  L\,  and  L2  be  Z)-dimensional  lattices  such  that 
L\  -<  L2  -<  Lq.  Consider  the  short  exact  sequence 

0  — >  B  — >  A  — ►  A/fl  — >  0, 

where  A  =  L0/Li  and  5  =  L2/Li  -<  A.  If  G  C  L0  is  a  complete  system  of  residues  for 
A,  Kp  C  L2a.  complete  system  of  residues  for  B,  and  GQ  C  £0  such  that  G  =  Kp  +  GQ 
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and  KpDGa  =  {0},  then  for  every  r  €  G  there  exist  unique  rp  €  Kp  and  ra  €  Ga  such 
that  r  =  rp  +  r0.  Furthermore,  card(G)  =  n,  card(A^j)  =  n^,  and  card(GQ)  =  na. 

Example.  Let  LQ  =  Z,L2  =  10Z,  and  Ij  =  1000Z.  Then  A  =  Z/1000Z  and 
5  =  10Z/1000Z.  The  set  G  =  {0, 1, ... ,  999}  C  Z  is  a  complete  system  of  residues 
for  A.  Furthermore,  if  Kp  =  {0, 10, 20, ... ,  990}  C  10Z  and  Ga  =  {0, 1, . . . ,  9},  then 
Kp  is  a  complete  system  of  residues  for  B,  Kp  C\Ga  =  {0}  and  G  —  Kp  +  Ga- 
in the  frequency  domain,  we  have  L*0  X  L\  <  L\  (Corollary  16).  Consider  the 
short  exact  sequence 

0  — ►  B'  — ►  A'  — y  A' IB'  —*  0, 

where  A'  =  L*JL'Q  and  5'  =  L"2/L'0  -<  A'.  Note  that  B'  £  L\j ' L\.  (To  understand  the 
relationship  between  B  and  B\  let  r  =  r  +  Li  €  B  =  L2/Li  and  5  =  s  +  LJ  €  B'  = 
L^/Lq.  Since  r  €  L2  and  s  €  L^  we  nave  (r> 5)  e  ^-  Hence,  (r,  5)  =  exp[-27n'(r,  s)]  = 
1.  Thus,  r£B  and  s  €  B'  implies  (r,s)  =  1.)  If  G'  C  IJ  is  a  complete  system  of 
residues  for  A',  K'a  C  L"2  a  complete  system  of  residues  for  S',  and  G^  C  L\  such 
that  A'a  n  G'p  =  {0}  and  G'  =  K'Q  +  G'p,  then  for  every  5  €  G'  there  exist  unique 
5Q  €  K'a  and  s^  €  G^  such  that  s  -  sa -\-  sp.  Corollary  16  and  Proposition  18  imply 
B'  =  A/B  and  A'/B'  S  B.  Since  A'  =  A,  we  have  card(G')  =  n,card(AQ)  =  na, 
and  card(G^)  =  np  =  c»xd(A'/Bt). 

Example.  If  L0  =  Z,L2  =  10Z,  and  L,  a  1000Z,  then  L'0  =  Z,L\  =  0.1Z, 
and  LJ  =  0.001Z.  Furthermore,  A'  =  0.001Z/Z  and  5'  =  0.1Z/Z.  The  set  G'  = 
{0,0.001,0.002,  ...,0.999}  C  0.001Z  is  a  complete  system  of  residues  for  A'  and 
K'a  =  {0,0.1,0.2,  ...,0.9}  C  0.1Z  is  a  complete  system  of  residues  for  B'.  If  G'p 
=  {0,0.001, 0.002,...,  0.099}  C  0.001Z,  then  K'a  n  G'p  =  {0}  and  G'  =  K'Q  +  G'p. 
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If  a  €  L(A)  and  a'  =  Ta{o.)  €  L(A'),  then  the  spatial  domain  CT  factoriza- 
tion for  quotient  lattices  is  determined  by 

a'(s)  =  -—   V  —=   Y]  a(r)(r0,s0)  I  ■  {ra,s0)  )  •  (ra,sa),  (3.6) 

for  all  s  =  sa  +  Bfi  e  K'Q  +  G'0  =  G",  where  r  =  r0  +  ra  6  K0  +  Ga  =  G.  If  a'  €  L(A') 
and  a  =  ^*71(o')  €  L(A),  then  the  inverse  or  frequency  domain  CT  factorization  for 
quotient  lattices  is  determined  by 


k&ifas,^ 


a{r)  -  —    2^  -^    }_   a'(s)(rQ,sQ)  )  •  {ra,s0)  J  •  (r^s/j), 

for  all  r  =  rp  +  ra  £  K0  +  Ga  =  G,  where  s  =  sa  +  s0  €  K'a  +  (?fi  =  G'. 


CHAPTER  4 
PERMUTOHEDRAL  AGGREGATES 


4.1     Introduction 


The  purpose  of  this  chapter  is  to  present  a  brief  introduction  to  certain  hier- 
archical arrangements  of  points  (or  cells)  known  as  permutohedral  aggregates.  These 
aggregates  and  their  geometric  duals  are  the  objects  of  interest  in  this  dissertation. 
This  chapter  applies  the  material  from  Chapter  3  to  quotient  lattices  that  are  as- 
sociated with  permutohedral  aggregates.  Specifically,  we  describe  the  Cooley-Tukey 
factorization  algorithms  in  the  spatial  and  frequency  domains,  respectively,  for  these 
quotients. 

Section  4.2  gives  a  geometric  description  of  permutohedral  aggregates.  In 
dimensions  two  and  three  they  are  hexagonal  and  truncated  octahedral  aggregates, 
respectively.  See  Conway  and  Sloane  [4]  for  the  necessary  background.  Section  4.3 
discusses  quotient  lattices  of  the  form  A\  =  L/NXL,  where  L  is  a  D-dimensional 
lattice  that  has  a  ring  structure,  N  is  an  endomorphism  of  L,  X  is  a  positive  integer, 
and  N  L  =  N  (L).  A  A-level  permutohedral  aggregate  is  a  particular  set  of  coset 
representatives  of  the  quotient  A\.  These  representatives  are  discussed  in  greater 
detail  in  Section  4.4.  Section  4.5  contains  statements  of  the  Cooley-Tukey  (CT) 
factorization  and  inverse  CT  factorization  in  terms  of  coset  representatives  of  A\  and 
the  geometric  duals  A'x.  Although  the  quotient  lattices  A\  have  a  ring  structure,  the 
factorization  algorithms  discussed  here  do  not  depend  on  the  multiplicative  structure 
of  these  quotients. 
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4.2     Geometric  Description 

Permutohedral  aggregates  can  be  defined  constructively  starting  with  a  poly- 
tope,  known  as  a  permutohedron,  associated  with  the  classical  Z)-dimensional  root 
lattice 

Ld  =  {(xo,xi,...,xD)  €  ZD+1  :  x0  -\ (-  xD  =  0}, 

which  lies  in  the  hyperplane  ^  xt  =  0  in  RD+X.  A  /^-dimensional  permutohedron  is 
the  Voronoi  cell  of  the  dual  lattice  /£,;  its  vertices  are  the  (D  +  1)!  points  obtained 
by  permuting  the  coordinates  of 


( 


\(-D),1-(-D  +  2),l-(-D  +  4),...,1-(D-2),1-D 


Permutohedra  tile  /^-dimensional  Euclidean  space. 

In  dimension  one,  L\  and  L\  can  be  represented  by  Z  in  R.  The  Voronoi 
cells  are  integer  translates  of  the  interval  [—1/2, 1/2].  A  level  1  aggregate  consists  of 
the  points  G  =  {-1,0, 1}.  For  A  €  Z+,  a  A-level  aggregate  is  the  set  of  3A  integers 
G\  —  YlkZo  %kG.  We  also  refer  to  the  set  of  Voronoi  cells  associated  with  G\  as  a 
A-level  aggregate. 

In  dimension  two,  Li  and  L\  are  hexagonal  lattices.  The  hexagonal  lattice  is 
an  optimal  packing  strategy  for  disks  in  the  plane.  The  Voronoi  cells  of  a  hexago- 
nal lattice  are  regular  hexagons.  A  level  1  hexagonal  aggregate  consists  of  a  central 
hexagon  together  with  six  neighboring  hexagonal  tiles  as  shown  in  Figure  4.1.  A 
level  2  hexagonal  aggregate  consists  of  seven  neighboring  level  1  aggregates.  In  this 
manner,  higher-level  aggregates  are  constructed  from  lower-level  aggregates.  The  hi- 
erarchical structure  of  a  level  3  hexagonal  aggregate  consisting  of  73  =  343  hexagonal 
cells  is  shown  in  Figure  4.2. 

In  dimension  three,  L3  is  a  face-centered  cubic  (fee)  lattice  and  L\  is  a  body- 
centered  cubic  (bec)  lattice.  The  fee  lattice  is  an  optimal  packing  strategy  for  spheres 
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Figure  4.1:  Level  1  hexagonal  aggregate. 
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Figure  4.2:  Hierarchical  structure  of  a  level  3  hexagonal  aggregate. 
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of  equal  radius  in  three-space.    The  Voronoi  cells  of  an  fee  lattice  are  rhombic  do- 
decahedra  (Figure  4.3).    The  Voronoi  cells  of  a  bec  lattice  are  truncated  octahedra 


Figure  4.3:  Rhombic  dodecahedron. 

(Figure  4.4).  The  bec  lattice  can  also  be  defined  as  the  set  of  vectors  (x,y,  z)  €  K3 
where  the  coordinates  x,  y,  and  z  are  all  even  or  all  odd  integers.  Figure  4.5  shows  a 
level  1  truncated  octahedral  aggregate,  which  consists  of  a  central  truncated  octahe- 
dron and  14  neighboring  cells.  Incidentally,  beehive  chambers  are  hexagonal  prisms 
capped  by  rhombic  dodecahedral  ends  (see  Kaproff  [12]). 

In  the  next  section  we  give  an  alternate  description  of  permutohedral  aggre- 
gates. 

4.3     Quotient  Lattices  Associated  with  Permutohedral  Aggregates 

Let  p  represent  the  polynomial  p(x)  =   1  +  x  +  •  ■  ■  +  xD,  where  D  €  Z+, 
and  let  A  =  Z[x]/(p),  a  commutative  ring  with  identity  T  =  1  +  (p).    Let  £  =  x  = 
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Figure  4.4:  Truncated  octahedron. 


x  -f  (p)  and  let  r  =  2  —  f .  The  quotient  ring  A  is  a  free  abelian  group  of  rank  D 
with  ordered  basis  X  =  (1,£,...  ,  fD-1),  where  1  is  the  identity  element.  We  have: 

££'  =  — (l  +  £  + y  £Z)_1)  and  (D+1  =  1.  The  quotient  ring  A  can  be  realized  as  a 

/^-dimensional  lattice  in  RD  by  mapping  the  basis  elements  to  any  set  of  D  linearly 
independent  vectors  in  RD.  We  will  simply  identify  the  basis  elements  in  X  with  a 
chosen  set  of  D  linearly  independent  vectors  in  WLD  and  regard  the  elements  of  X  as 
labels  or  names  for  the  assigned  vectors.  The  element  r  determines  an  endomorphism 
NT  on  A  defined  by  NT  :  p  i->  rp,  for  all  p  €  A.  For  example,  if  D  =  2,  then 
A/T(£)  =  T£  =  (2- £){  =  2£ - ?  *  2£  + 1  +  {  =  1  +  3f  Let  NTA  =  NT(A).  Note  that 
NTA  m  tA  =  (r)  is  a  principal  ideal  of  A.  Furthermore,  tX  =  (r,  r£, . . . ,  t^-1)  is  a 
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Figure  4.5:  Level  1  truncated  octahedral  aggregate. 
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basis  of  tA.  With  respect  to  the  basis  X, 

/     2 

-1 


NT  = 


0    0 
2    0 
0-12 


0 


0    0 

0    0 


...     0 

1  \ 

...     0 

1 

...      0 

1 

•••       2 

1 

...    _i 

3/ 

(4.1) 


Furthermore,  det  JVT  =  2D+1  -  1  (see  Vince  [22]).  More  generally,  for  A  €  Z+, 
NXA  =  A^(A)  =»  txA  =  (rA),  is  an  ideal  of  A,  and  det  Nx  =  nx.  Kitto,  Vince,  and 
Wilson  [13]  prove  the  following  isomorphism  of  rings. 

Theorem  19  If  D,\  6  Z+,  then  A/NXA  =  Zn\  (ring  isomorphism),  where  n  = 
2D+1  -I.  In  particular,  A/NXA  is  cyclic  and  card(A/NxA)  =  nx. 

Let  {uj  :  j  =  0, 1, . . . ,  D  —  1}  be  a  set  of  D  linearly  independent  vectors  in  RD 
and  denote  by  U  the  ordered  basis  (uo,  u\, . . . ,  uq-i)  of  Rc.  Identify  £J  with  Uj,  for 
j '  =  0, 1, . . . ,  D  —  1,  and  denote  by  L  the  jD-dimensional  lattice  generated  by  U.  Let 
N  be  the  endomorphism  of  L  such  that  the  matrix  of  Af  with  respect  to  the  basis  U  is 
given  by  the  matrix  in  equation  4.1.  Thus,  AX  =  A^(L)  is  a  /^-dimensional  sublattice 
of  L  and  the  index  of  NL  in  L  is  n  =  det  A^  =  2D+l  -  1.  More  generally,  if  A  €  Z+ 
and  NXL  =  NX(L),  then  NXL  is  a  .D-dimensional  sublattice  of  L  that  is  generated 
by  the  set  of  vectors  {Nx(uj)  :  j  =  0, 1, . . . ,  D  —  1}.  Furthermore,  the  index  of  NXL 
in  L  is  det  A^  =  (det  N)x  =  nx. 

We  turn  our  attention  to  quotient  lattices  derived  from  L  and  the  sublattices 
NXL.  Let  A  =  L/NL,  that  is,  A  is  the  cokernel  of  N  e  End(L),  the  set  of  endo- 
morphisms  of  L.  We  have  card(/l)  =  n  —  2D+i  —  1.  More  generally,  for  A  6  Z+,  let 
Ax  =  L/NXL,  that  is,  A\  is  the  cokernel  of  Nx  €  End(L).  We  have  card(AA)  =  nA. 

Remark  2  A  X-level  permutohedral  aggregate  (in  the  spatial  domain)  is  a  particular 
set  of  co set  representatives  of  A\  as  a  quotient  lattice.  For  permutohedral  aggregates, 
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the  Uj  €  U  are  chosen  such  that  the  set  of  D  +  1  points  {uj  :  j  =  0, 1, . . . ,  D},  where 
U£>  =  —  Y^°=o  ui>  ls  ^e  set  of  vertices  of  a  regular  D-simplex  with  barycenter  at  the 
origin.  See  Lucas  [14]  for  a  particular  choice  of  coordinates  for  these  vectors. 

By  the  lattice  duality  relationship  (Proposition  15),  NL  -<  L  implies  L*  -< 
(NL)* .  Hence,  the  geometric  counterpart  of  A  in  the  frequency  domain  is  the  quotient 
lattice  A'  =  (NL)*/L*.  Similarly,  for  A  6  Z+,  the  geometric  counterpart  of  A\  in  the 
frequency  domain  is  A'x  =  (NXL)*/L'.  We  can  regard  the  (injective)  endomorphisms 
Nx  as  .//-invariant,  invertible,  linear  operators  on  RD  (and  therefore  bounded).  Recall 
that  every  bounded  linear  operator  on  a  Hilbert  space  has  a  unique  adjoint. 

Theorem  20  If  H  is  a  Hilbert  space  and  M  is  a  bounded  linear  operator  on  H, 
then  there  exists  a  unique  linear  operator  M." ,  known  as  the  adjoint  of  M.,  satisfying 
{M(x),y)  =  (x,M'(y)),forallx,ye  H. 

For  the  elementary  results  from  linear  operator  theory  in  this  section  and  for  a 
proof  of  Theorem  20,  which  uses  the  Riesz  representation  theorem  for  a  Hilbert  space, 
see  DeVito  [6].  We  shall  use  the  following  results  for  a  linear  operator  M.:  if  M.  is 
an  invertible  linear  operator  and  M*  is  the  adjoint  of  M ,  then  M*  is  invertible  and 
{M*)~l  =  («M-1)*;  if  M  is  given  by  a  real-valued  square  matrix  with  respect  to  a 
particular  basis,  then  M*  =  M\  the  transpose  matrix  with  respect  to  the  same  basis; 
if  A  e  Z+,  then  {Mx)m  =  (M*)x;  finally,  detAT  =  detM.  The  next  proposition  will 
allow  us  to  express  A'  and  A'x  in  terms  of  W,  the  unique  adjoint  of  the  L-invariant, 
linear  operator  ./V.  In  general,  for  /  a  .D-dimensional  lattice  and  AA  an  /-invariant 
linear  operator  on  RD,  let  MI  =  M(I).  Note  that  if  M  is  invertible,  then  MI  is  a 
Z)-dimensional  lattice  and  MI  -<  I  (sublattice  and  subgroup  relationship). 
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Proposition  21  ///  is  a  D -dimensional  lattice,  I'  is  the  dual  lattice  of  I,  M  is  an 
I -invariant,  invertible  linear  operator  on  RD,  and  M"  is  the  unique  adjoint  of  M, 
then  (1)  (MI)'  =  (M')~lI',  and  (2)  if  V  =  {MI)',  then  M"  is  I' -invariant. 

Proof.  Part  (1).  We  first  show  (MI)'  C  (A**)"1/*.  Let  u  €  (MI)'.  Then 
(M(r),u)  =  (r,M'(u))  €  Z,  for  all  r£/.  Hence,  M'(u)  €  /*,  that  is,  there  exists 
s  6  /'  such  that  M'(u)  =  s,  or  u  =  {Mm)"l(s).  Thus,  u  €  (M')~lI'  and  (MI)'  C 
(M*)~lI'.  To  show  (M*)~lI'  C  (MI)*,  let  s  €  /*•  We  must  show  (M*)~l(s)  € 
(MI)*,  that  is,  ^(^(M*)"^))  €  Z,  for  all  r  €  /.  Let  r  €  /.  Then  X(r)  G  X/ 
and  {M{r),(Mm)'l{*))  =  (^M^-1)*!5))  =  (M'^MW)**)  =  (r,5)  6  Z,  since 
r  €  /  and  s  6  /*.  Hence,  (AT)"1^)  6  (-^/)"  and  (M')~lI'  C  (A<J)'. 

Part  (2).  The  hypotheses  imply  MI  <  I.  Let  /'  =  (MI)'  =  (AT)"1/*, 
using  the  first  part  of  the  proposition.  Then  7"  =  M'l'.  Now,  MI  -<  /  implies 
/•  -<  (MI)"  (Proposition  15),  that  is,  MT  <  V.  Hence,  M'  €  End(/').  □ 

For  D,\  €  Z+  fixed,  let  L\  =  (NXL)'.  Part  (1)  of  the  proposition  implies 
L'x  =  ((iVA)*)-1L*  =  (N')~XL';  part  (2)  implies  (Nx)*  =  (N')x  is  LA-invariant. 
Hence,  L'  =  (N')XL'X  and  A'x  =  L'X/(N')XL'X.  Just  as  Ax  is  the  cokernel  of  Nx  € 
End(L),  A'x  is  the  cokernel  of  (N*)x  €  End(Z/A).  The  corresponding  short  exact 
sequences — for  A\  and  A'x,  respectively — are 

0  — ►  NXL  — ►  L  — *  L/ATAL  — ►  0, 

in  the  spatial  domain,  and 

0  -4  (N')XL'X  -*  Z'A  — ►  L'X/(N')XL'X  -4  0, 

in  the  frequency  domain. 
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Remark  3  Note  that  N"  is  L'x-invariant.   To  see  this,  consider  the  chain  of  lattices 

NXL  <  Nx~lL  < <  NL  <  L.  In  particular,  NXL  -<  Nx~lL  implies  (A^A"1L)*  * 

(NXL)*,  or  (N*y-XL*  <  (N')-xLm,  by  the  first  part  of  Proposition  21.    Hence,  if 
L'x  =  (N')-XL',  then  N"L'X  -<  L'x. 

Denote  by  L(A\)  the  nA-dimensional  space  of  complex-valued  functions  on 
A\,  and  denote  by  L(A'X)  the  nA-dimensional  space  of  complex-valued  functions  on 
A'x.  Furthermore,  denote  the  discrete  Fourier  transform  (equation  3.1)  of  a  function 
a  €  L(A\)  by  a'  =  F\(a)  G  L(A'X),  and  denote  the  inverse  DFT  (equation  3.2)  of  a 
function  a'  e  L{A'X)  by  o  =  F^l{a')  €  L{A\). 

4.4     Aggregates  as  Sets  of  Coset  Representatives 

Let  /  be  a  Z)-dimensional  lattice  in  RD  and  M  €  End(7)  such  that  |  det./V{|  > 
1.  A  set  R  C  I  is  called  a  residue  system  for  M.  if  0  (E  R  and  /  =  (J{r  +  MI  :  r  €  R} 
(see  Bandt  [3]).  Thus,  R  is  a  set  of  coset  representatives  of  I/MI  and  card(i?)  = 
det./V4.  A  finite  address  of  a  lattice  point  r  £  /  is  a  finite  sequence  r\_i  •  •  -Tiro 
(A  €  Z+)  such  that  r*  £  R,  for  all  ifc,  and  r  =  £^j  A^fc(rit).  Vince  [22]  shows  that 
the  finite  address  of  a  lattice  point,  if  it  exists,  is  unique.  A  X-level  aggregate  (not 
necessarily  a  permutohedral  aggregate)  is  the  set  of  lattice  points  having  addresses  of 
length  no  greater  than  A.  We  also  refer  to  the  corresponding  set  of  Voronoi  cells  as  a 
A-level  aggregate.  The  following  lemma  states  that  A-level  aggregates  are  particular 
sets  of  coset  representatives  of  I / Mxl . 

Lemma  22  Let  D,  A  €  Z+.  7/7  is  a  D-dimensional  lattice,  M  G  End{I),  [  det  Af  |  > 
1,  and  R  is  a  residue  system  for  M}  then 

A-l 

Rx  =  J2Mk(R) 

k=0 
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is  a  set  of  coset  representatives  of  I /MI. 

Proof.     Since  R  is  a  set  of  coset  representatives  of  I /MI  and  MI  C\  R  =  {0},  we 
have  I  =  MI  +  R.  Now, 

/    =    MI+R 

=    M(MI+R)  +  R 
=    M2I  +  M(R)  +  R 

■ 

=    MXI  +  MX-1{R)  +  ---  +  M{R)  +  R 
=    MXI  +  RX. 

Hence,  R\  is  a  set  of  coset  representatives  of  I/MXI.  □ 

For  D  <E  Z+,  let  n  =  2D+1  -  1,  p(x)  =  1  +  x  +  •  •  •  +  xD,  A  =  Z[x]/(p),  and 
r  =  2  —  £,  where  £  =  x  +  (p),  as  before.  Kitto,  Vince,  and  Wilson  [13]  prove  the  next 
result  (we  omit  the  use  of  overbars  for  elements  of  A/NTA;  NT  is  the  endomorphism 
of  A  defined  in  Section  4.3). 

Lemma  23  Each  of  the  following  subsets  of  A  is  a  set  of  coset  representatives  of  the 
quotient  ring  A/NTA: 

1.  R  =  {0,l,...,n-1}; 

2.  S  =  {bD(D  +  ---  +  b1Z  +  b0:bJ  G  {0,1},  forj  =  0,l,...,D,    not  all  6;  =  1}. 

Note  that  elements  in  the  second  set  in  Lemma  23  can  be  written  as  strings  (of 

length  D+l)  of  binary  digits,  or  bits.  Thus,  we  can  represent  s  =  bD£D-\ h&if+&o  € 

S  by  the  string  (bD  ■  ■  ■  bib0).   Conversely,  if  (bD  ■  ■  ■  bib0)  is  the  binary  representation 
ofreR,  then  r  =  s  mod  NTA,  where  s  =  bD£D  +  ■  •  •  +  b^  +  b0  =  (bD  ■  ■  ■  bibo)  €  S. 


43 

Example.  If  D  =  1,  then  A  =  Z[x]/(1  +  x)  is  isomorphic  to  Z,  R  =  {0, 1,2}  and 
S  =  {0, 1,  £}.  Note  that  2  ee  £  mod  ATA. 

Example.  If  D  «  2,  then  A  =  Z[x]/(1  +  x  +  x2),  A  ■  {0,1,2,3,4,5,6}  and 
S  =  {0, 1,  f ,  £  + 1,  £2,  f2  +  £}.  The  following  equivalences  (modulo  A7,- A)  hold  between 
the  elements  of  R  and  the  elements  of  5: 

2  =  e, 

3    =    £  +  1, 

4  =  ea, 

5  =  e+h 

6  =  e+t. 

Note  that  the  binary  representation  of  6  6  R  is  (110),  which  represents  the  equiva- 
lence class  l-£2-|-l-£  +  0-l.  Similarly  for  the  other  elements  of  R. 

We  now  consider  sets  of  coset  representatives  of  the  quotients  A\  and  A'x.  Let 
U  =  (u0,  ui, . . . ,  ud-i)  be  an  ordered  set  of  D  linearly  independent  vectors  in  RD  and 
let  ud  =  —  5^7=0  Uj.  Under  the  identifications  £J  f->  Uj,  for  j  =  0, 1, . . . ,  D,  the  set 


<  Y^  biu3  •  h  €  {0, 1},  for  j  =  0, 1, . . . ,  D,  not  all  bj  =  1  >  (4. 


2) 


is  a  residue  system  for  N  :  L  — >  L,  where  L  is  the  /^-dimensional  lattice  generated 
by  £/  and  A  is  the  endomorphism  of  L  specified  by  the  matrix  in  equation  4.1  with 
respect  to  the  basis  U.  The  set  G  is  equivalent  to  the  second  set  in  Lemma  23  under 
the  identifications  &  «-*  Uj.  Recall  that  det  N  =  2D+1  -  1.  By  Lemma  22,  the  set 

A-l 


GX  =  £  Nk(G) 


k=0 
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is  a  set  of  coset  representatives  of  A\  =  L/NXL. 

In  the  frequency  domain,  let  U'x  =  (u'0,  u[, . . . ,  u'D_x)  be  an  ordered  set  of  D 
linearly  independent  vectors  in  RD  that  generate  the  Z)-dimensional  lattice  L'x  = 
(NxL)m  =  {N')-xLm,  and  let  u'D  =  -  JjJJ  uj  such  that 

G'=l^T  bju'j  :  b3  €  {0, 1},  for  j  =  0, 1, ... ,D,  not  all  6j :  =  1  I  (4.3) 

is  a  residue  system  for  N*  :  L'x  -+  L'x.  Note  that  det  N*  =  det  N.  By  Lemma  22,  the 

set 

A-l 


G'x  =  J2(N')k(G') 

is  a  set  of  coset  representatives  of  A'x  =  L'X/(N*)XL'X. 

Example.  Let  D  —  2  and  let  L  be  the  2-dimensional  lattice  in  the  spatial  domain 
generated  by£/  =  (u0,Ui)  and  let  u2  —  —(u0  +  Ui).  We  have  G  =  {62«2  +  &i«i  +&0W0  : 
ij  €  {0,1},  not  all  6,-  =  1}, 

with  respect  to  the  ordered  basis  U,  and  det  N  —  7  =  card(G).  If  A  =  2,  then  L'2 
is  the  lattice  (N2L)*  =  (N*)~2L'  in  the  frequency  domain.  Suppose  Z/2  is  generated 
by  U'2  =  KX)  and  let  u'2  =  -(u'0  +  u[).  Then  G'  =  {b2u'2  +  6^;  +  b0u'0  :  bj  € 
{0,1},  not  all  bj  =  1}  is  a  set  of  coset  representatives  of  A'2  =  L'2/(N*)2L'2.  If 
u0  =  (1,0)'  and  ui  =  (-1/2,  y/%f2)*,  then  I  is  a  hexagonal  lattice  in  the  plane 
(congruent  to  the  lattice  L2  in  Section  4.2);  the  associated  Voronoi  cells  are  hexagons. 
Figure  4.6  depicts  a  level  1  hexagonal  aggregate  (the  set  G)  in  the  spatial  domain. 
The  figure  shows  the  seven  lattice  points  and  the  corresponding  hexagonal  tiles.  Note 
that  the  ordered  pair  (u0,«i)  defines  a  right-handed  coordinate  system  in  the  plane. 
Level  two  hexagonal  aggregates  in  the  spatial  and  frequency  domains  are  shown  in 
Figures  4.7  and  4.8,  respectively.  Geometrically,  the  two  domains  differ  in  scale, 
orientation,  and  chirality  (or  handedness). 


1.5 


-1.5  -1  -0. 


-1.5 
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Figure  4.6:  Level  1  hexagonal  aggregate  in  the  spatial  domain. 


46 


Figure  4.7:  Level  2  hexagonal  aggregate  in  the  spatial  domain. 
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Figure  4.8:  Level  2  hexagonal  aggregate  in  the  frequency  domain. 
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Example.  If  D  =  3,  then 

/     2  0  1 

A=       -1  2  0 

\     0  -1  3 

and  detTV  =  15  =  card(G).  Choose  u0  =  (1,0,0)',U!  =  (-l/3,2v/2/3,0)',u2  ■ 
(— 1/3,  — \/2/3,2/\/6)'.  These  vectors  generate  a  body-centered  cubic  lattice  (con- 
gruent to  the  lattice  L^  from  Section  4.2).  The  associated  Voronoi  cells  are  truncated 
octahedra.  A  level  1  truncated  octahedral  aggregate  is  shown  in  Figure  4.5. 


4.5     The  Cooley-Tukey  Factorization  for  Permutohedral  Aggregates 

Fix  D,\eZ+  and  let  n  =  2D+l  -  1. 

In  the  spatial  domain,  consider  the  short  exact  sequence 

0  — ►  Bp  — ►  Ax  — ►  Qa  — ►  0, 

where  Ax  =  L/NXL,  Bp  is  the  subgroup  NaL/NxL  of  Ax,  and  £„  =  Ax/Bp.  For 
1  <  7  <  A,  let  G-,  =  YllZo  Nk(G),  where  G  is  a  residue  system  for  the  endomorphism 
N.L-+L.ttKfi=  Na(G0),  then  Kp  nGa  =  {0},  and  GX  =  KP  +  Ga  is  a  complete 
set  of  coset  representatives  of  A\.  Hence,  for  every  r  £  Gx,  there  exist  unique 
rp  E  Kfi  and  ra  €  GQ  such  that  r  =  rp  +  ra.  We  regard  Kp  as  the  set  of  high-order 
representatives  in  Ga;  similarly,  we  regard  Ga  as  the  set  of  low-order  representatives 
inGA. 

In  the  frequency  domain,  let  L'x  =  (NXL)*  and  consider  the  short  exact  se- 
quence 

where  A'x  =  L'X/(N*)XL'X,  B'Q  is  the  subgroup  {N'fL'J{N')xL'x  of  A'x,  and  <%  = 
A'JB'a.  For  1  <  7  <  A,  let  G\  =  £2:J(A*)*(G"),  where  G"  is  a  residue  system  for 
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the  endomorphism  N'  :  L'x  -¥  L'x.  If  K'a  =  (N'f(G'a),  then  K'anG'0  =  {0},  and 
G'x  =  K'a  +  G'p  is  a  complete  set  of  coset  representatives  of  A'x.  Hence,  for  every 
s  €  G'x,  there  exist  unique  sa  €  K'a  and  sp  €  G0  such  that  s  =  sa  +  5/?-  We  regard 
A'q  as  the  set  of  high-order  representatives  in  G'x,  and  G'0  as  the  set  of  low-order 
representatives  in  G'x. 

If  a  £  L(AX),  then  the  discrete  Fourier  transform  a'  =  T\{a)  €  L{a)  can 
be  computed  by  the  spatial  domain  Cooley-Tukey  factorization  (Theorem  13  and 
equation  3.6)  as 

a'^  =  ^  E    b^E  «(»•)(**.«/>)]  ■('«*))  •(»•«,««), 

for  all  s  =  sa  +  s0  e  K'a  +  G'p  =  G'x,  where  r  =  r0  +  ra  e  K0  +  Ga  =  Gx.  If 
a'  6  L(A'X),  then  the  inverse  DFT  a  =  ^l(o')  €  I(Aa)  can  be  computed  by  the 
frequency  domain  CT  factorization  (Corollary  14  and  equation  3.7)  as 


J^    a'(s)(ra,sa)      ■  (ra,sp)      ■  (rp,sp), 


for  all  r  =  r0  +  raeK0  +  Ga  =  G\,  where  s  =  sa  +  S0  6  K'a  +  G'0  =  G'x. 

The  ideals  B0  do  not  split  in  A\.  For  example,  in  dimension  two,  A\  =  Z/7AZ, 
which  does  not  split  since  7  is  prime.  Thus,  the  Good-Thomas  algorithm  (as  described 
by  An  et  al.  [1],  for  example)  is  not  applicable  to  discrete  functions  defined  on  the 
quotients  .4^-  The  Good-Thomas  algorithm  is  related  to  the  Chinese  remainder 
theorem  and  is  used  to  compute  the  DFT  of  functions  with  rectangular  domains,  that 
is,  domains  of  the  form  Z„,  x  •  •  •  x  ZnD.  The  algorithm  exploits  the  linear  separability 
of  the  Fourier  transform  by  first  computing  the  one-dimensional  transform  over  the 
rows  of  an  image  and  then  computing  the  one-dimensional  transform  over  the  columns 
of  the  resulting  intermediate  image  (similarly  for  higher  dimensions). 


CHAPTER  5 
THE  GENERALIZED  MATRIX  PRODUCT 


5.1     Introduction 

The  goal  of  this  chapter  is  to  introduce  the  generalized  matrix  product  over 
digital  indexing  sets.  The  generalized  matrix  product,  formulated  by  Ritter  [17], 
generalizes  standard  matrix  multiplication.  By  using  this  product,  data  reordering 
stages  in  the  FFT  algorithm  are  eliminated.  Data  reordering  is  a  feature  of  stan- 
dard FFT  algorithms.  Zhu  and  Ritter  [26]  examined  the  computational  advantages 
of  using  the  generalized  matrix  product  in  signal  processing  applications.  Digital 
indices  address  the  elements  of  permutohedral  aggregates  and  their  geometric  duals. 
This  form  of  addressing  is  known  as  generalized  balanced  ternary  (GBT)  and  was 
developed  by  Gibson  and  Lucas  [8].  The  generalized  matrix  product  described  by 
Ritter  and  Zhu  is  not  easily  adapted  to  data  arrays  indexed  by  strings  of  digits,  as  in 
the  GBT  scheme  for  permutohedral  aggregates.  While  the  generalized  matrix  prod- 
uct and  GBT  addressing  for  permutohedral  aggregates  have  been  studied  previously 
and  independently,  the  product  defined  in  this  chapter  is  tailored  to  processing  data 
arrays  that  are  associated  with  permutohedral  aggregates  and  their  duals. 

Let  ¥  be  a  ring  and  denote  by  F;Xm  the  set  of  /  x  m  matrices  with  entries  in  ¥. 
The  definition  of  standard  matrix  multiplication,  which  is  a  map  F/Xm  xFmxr  — »  F/Xr, 
requires  the  number  of  columns  in  the  left  factor  to  equal  the  number  of  rows  in  the 
right  factor.  With  p  €  Z+,  the  generalized  matrix  product  of  order  p  (or  p-product) 
is  a  map  ©p  :  F,xmp  x  Fp,xr  -4  Flqxmr.  Hence,  the  p-product  A  0P  B  requires  that 
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p  divide  the  number  of  columns  in  the  left  factor  and  the  number  of  rows  in  the 
right  factor.  Note  the  dimensions  of  the  resulting  matrix.  We  will  also  make  use 
of  the  Hadamard  product,  a  componentwise  operation  on  matrices  having  the  same 
dimensions,  defined  as  follows.  For  A,B  e  F;xm,  the  Hadamard  product  C  =  (efc)  of 
A  —  (a,j)  and  B  =  (6,j),  denoted  by  A  o  B,  is  given  by  C{j  =  aijbij,  for  all  indices  i 
and  j.  Hence,  the  Hadamard  product  is  a  map  o  :  F/Xm  x  F/Xm  — ►  F/Xm. 

Section  5.2  defines  T -indices,  indexing  sets  consisting  of  finite  strings  of  digits 
that  address  the  elements  of  permutohedral  aggregates  and  their  geometric  duals 
(that  is,  coset  representatives  of  A\  and  A'x,  respectively).  Addressing  with  T-indices 
can  be  implemented  in  terms  of  strings  of  bits,  or  binary  digits,  and  yields  very 
efficient  data  processing  algorithms  (see  Gibson  and  Lucas  [8],  and  Kitto,  Vince,  and 
Wilson  [13]).  Section  5.3  defines  the  generalized  matrix  product,  or  p-product,  over 
T-indices,  and  Section  5.4  covers  digital  addressing  for  permutohedral  aggregates  and 
their  duals. 

5.2     Digital  Indexing  Sets 

The  elements  of  .D-dimensional  permutohedral  aggregates  and  their  duals  can 
be  addressed  using  finite  strings  of  digits.  We  refer  to  the  underlying  digital  indexing 
sets  as  T-indices,  defined  as  follows. 

Definition  4  (Gamma  Indices)  Let  D  €  Z+  and  n  =  2D+1  -  1. 

1.  Define  the  set 

r1  =  {0,l,...,n-l} 

together  with  the  order  relationship 

0-<l-<--<n-l; 
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2.  For  A  (E  Z+,  define  the  set 

Tx  =  {gx-i9\-2---9i9o  -9k  €  Ti,  for  k  =  0, 1, . . . ,  A  -  1} 

and  extend  the  order  relationship  on  T\  to  lexicographic  order  on  T\,  where  the 
rightmost  digit  g0  is  the  low-order,  that  is,  least  significant,  digit; 

3.  For  A, a,/?  €  Z+  such  that  A  =  a  +  j3,  define  the  set 

rar/3  =  {gh:gerQ,i3er0}. 

We  refer  to  the  sets  T\  as  Y-indices. 

Remark  4   The  following  properties  of  T-indices  follow  directly  from  Definition  4- 
Fix  D  €  Z+  and  let  n  =  2D+1  -  1. 

1.  The  set  T\  has  nx  elements; 

2.  If  A,  q,  /3  e  Z+  such  that  A  =  a  +  /?,  then 

TaTp  =  T\  =  FpFa. 


5.3     The  Generalized  Matrix  Product  for  Permutohedral  Aggregates 

The  generalized  matrix  product,  or  p-product,  over  T-indices  has  a  particularly 
simple  form.  Fix  D  e  Z+  and  let  n  =  2D+1  -  1.  Let  F  be  any  ring  and  a,0  €  Z+. 
Denote  by  FpQ  xr>  the  set  of  all  na  xn0,  F- valued  matrices  with  rows  indexed  by  the 
set  of  indices  Ta,  and  columns  indexed  by  the  set  of  indices  Tp.  In  this  notation, 
*  €  FrQxi  is  a  column  vector  of  length  na  indexed  by  rQ,  and  y  €  Fixr„  is  a  row 
vector  of  length  n0  indexed  by  T0.  If  transposition  is  denoted  by  superscript  t,  then 
x  €  FrQxi  implies  x*  eFlxra. 


53 

Definition  5  (Generalized  Matrix  Product  Over  Gamma  Indices)  Let  ¥  be 

any  ring,  D  G  Z+,  and  n  =  2D+1  -  1.  Let  A  €  FrQxr>  where  a,  /?  €  Z+,  and 
B  €  Fr^xTi,  where  7,0"  G  Z+.  //A  €  Z+  suc/i  Ma£  A  <  {/3,t},  Men  Me  p-product  of 
A  and  B,  with  p  =  nx,  is  the  na+1~x  x  n0~x+s  matrix 

C  =  ABpB  €  Fr.r^xr^r, 

defined  by 

C{gt,rm)  =  ^  ,%,rs)£(^,m), 
«erA 
/or  all  gt  €  ror7_A,  rm  €  IV  a  IV 

Remark  5  For  /u//  matrices  A  and  B  as  in  Definition  5,  computation  of  the  nx- 
product  C  =  A  ©„a  B  requires  na+0+-r+&-i  multiplications.  The  case  /?  =  A  =  7 
corresponds  to  standard  matrix  multiplication. 

Example.  If  D  =  1,  then  n  =  3  and  I\  =  {0, 1,2}.  If  x  €  Rlxr3  and  y  <E  Rr2xi, 
then 

x  =  (:r(000),  x(001),  x(002),  x(010), . . . ,  i(221),  x(222)) 
and 

y  =  (2/(00),  2/(01),  y(02),  2/(10), . . . ,  2/(21),  2/(22))'. 
If  A  =  2  and  z  =  x  ®nx  y  =  x  ©32  y,  then  2  €  lixr,  and 

z(r)  =  J2  x(™)y{s), 

s€T2 

for  r  =  0,1,2. 

Example.  If  L>  =  2,  then  n  =  7  and  I\  =  {0, 1, ...  ,6}.  Suppose  A  €  Cr2Xr3  and 
£eCr4xrv  IfC  =  A®7B(\  =  l),  thenCeCr5Xr2  and 


C{gt,rm)  =  J^  A(g,rs)B(st,m), 


*er, 
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for  all  (gt,rm)  €  r2r3  x  I^Tx  =  T5  x  T3. 

5.4     Digital  Addressing  for  Permutohedral  Aggregates 

Our  goal  in  this  section  is  to  express  the  elements  of  a  permutohedral  aggregate 
(and  of  the  corresponding  geometric  dual)  entirely  in  terms  of  digits,  or  T-indices.  We 
will  need  to  distinguish  between  digital  addresses  for  elements  of  permutohedral  ag- 
gregates in  the  spatial  domain  and  digital  addresses  for  elements  of  the  corresponding 
duals  in  the  frequency  domain.  To  this  end,  let  T[  =  Tj  and  F'x  =  T\  (A  6  Z+). 

Fix  D,  A  €  Z+.  Let  L  be  the  Z)-dimensional  lattice  generated  by  U  = 
(uq,  Ui,.  . .  ,ud-i)  and  let  N  be  the  endomorphism  of  L  specified  by  the  D  x  D 
matrix  in  equation  4.1  with  respect  to  the  basis  U.  Let  li'x  =  (u'0,  u[, . . . ,  U/j-i)  be  an 
ordered  set  of  D  linearly  independent  vectors  that  generate  the  lattice  L'x  =  (N*)~  L. 
Let  uq  =  —  X^,=o  uj  and  u'd  ~  ~  J2j=o  u'j-  Recall  the  sets  of  coset  representatives 
G,G',G\,  and  G'x  described  in  Section  4.4.  Note  that  any  element  of  Ti,  that  is, 
any  digit,  can  be  expressed  in  the  form  bo2D  +  •  •  •  +  bx2  +  bo,  where  the  bj  €  {0, 1}, 
for  all  j,  and  not  all  b3   =   1.     Furthermore,  the  bj  are  unique.    Assign  the  digit 

bn2D  H \-bi2  +  b0  6  Ti  to  the  vector  bouo-\ hfeiWi  +  60uo  €  G.  This  assignment 

of  digits  establishes  a  one-to-one  correspondence  between  Ti  and  G.  Now,  assign  the 

string  of  digits  gx-i  •  •  ■  gig0  €  Tx  to  the  vector  r  =  Nx~1(rx_1)  + hA^(r1)  +  r0  €  Gx, 

where  <ft  6  Ti  is  the  digit  corresponding  to  the  vector  r*  €  G,  for  k  =  0, 1, . . . ,  A  —  1. 
In  this  manner,  we  establish  a  one-to-one  correspondence  between  the  digital  in- 
dexing set  Fx  and  G\  (a  permutohedral  aggregate  in  the  spatial  domain).  Write 
9k  =  Tk  whenever  <ft  €  I\  is  the  digit  corresponding  to  r*  €  G  and  refer  to  the  string 
9\-i  ■  -gigo  =  Nx~1(gX-i)  H 1-  N(gi)  +  g0  as  the  digital  address  of  r. 

In  the  frequency  domain,  assign  the  digit  bD2D  -\ (-  b^  +  bo  G  T[  to  bDuo  + 

h  biUi  +  b0u0  €  G'.  This  assignment  of  digits  defines  a  one-to-one  correspondence 
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between  Ti  and  G' .  Assign  the  string  of  digits  /ia-i  ■  ■  ■  hxhQ  6  V'x  to  the  vector 
s  =  Nx~1(sx-i)  +  •  •  •  +  N{s\)  +  s0  €  G^,  where  /i;t  €  Tj  is  the  digit  assigned  to 
Sk  €  G',  for  k  =  0,1,...,  A  —  1.  This  assignment  of  digital  strings  defines  a  one- 
to-one  correspondence  between  the  digital  indexing  set  T'A  and  G'x  (a  permutohedral 
aggregate  in  the  frequency  domain).  Write  hk  =  sk  whenever  hk  €  T[  is  the  digit 
assigned  to  sk  €  G';  the  string  /ia-i  •  •  •  hih0  =  ^2kZo(N*)k(hk)  is  the  digital  address 
of  s. 

Recall  from  Section  4.4  that  the  finite  address  of  a  lattice  point,  if  it  exists, 
is  unique.  Hence,  the  digital  addressing  scheme  we  have  described  assigns  unique 
addresses  to  the  elements  of  permutohedral  aggregates  in  the  spatial  domain;  similarly 
for  their  counterparts  in  the  frequency  domain.  Also,  note  that  the  elements  of  V\ — 
or  the  corresponding  vectors  in  G — can  be  expressed  as  bit  patterns  of  the  form 
(bp  ■  •  •  6i60).  Similarly  for  the  elements  of  T[  and  the  corresponding  vectors  in  G'. 

Example.  Let  D  =  2  and  A  =  3.  The  digital  address  of  the  vector  N2(u2  +  u0)  + 
N(u2)  +  u0  =  A72(l-u2  +  0-tx1  +  l-uo)  +  A;(l-u2-r0-u1-r0-Uo)  +  (0-u2-|-0-u1-|-l-Uo)  = 
N2(101)  +  Ar(100)  +  (001)  =  N2(5)  +  N{4)  + 1  G  G3  is  541  €  [».  Note,  however,  that 
541  e  H,  is  the  digital  address  of  the  vector  (N')2{u'2  +  u'0)  +  N'(u2)  +  u'0  6  G'3  in 
the  frequency  domain. 

Example.  Figure  5.4  shows  an  indexed  level  3  hexagonal  aggregate  in  the  spatial 
domain. 
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Figure  5.1:  Indexed  level  3  hexagonal  aggregate  in  the  spatial  domain. 


CHAPTER  6 
DIGITAL  ALGORITHMS  FOR  PERMUTOHEDRAL  AGGREGATES 


6.1     Introduction 

This  chapter  brings  together  the  concepts  developed  in  earlier  chapters  to 
give  statements  of  the  discrete  Fourier  transform,  the  Cooley-Tukey  factorization, 
and  the  fast  Fourier  transform — together  with  the  respective  inverse  transforms — for 
complex-valued  functions  defined  on  permutohedral  aggregates  and  their  geometric 
duals.  The  algorithms  in  this  chapter  are  expressed  in  terms  of  the  generalized  matrix 
product  over  digital  indexing  sets. 

The  two  theorems  in  Section  6.2  define  the  discrete  Fourier  transform  (DFT) 
and  the  inverse  DFT  for  permutohedral  aggregates.  The  theorems  and  lemmas  in 
Section  6.3  are  digital  Cooley-Tukey  (CT)  factorization  algorithms  for  permutohedral 
aggregates.  Section  6.4  gives  statements  of  the  fast  Fourier  transform  (FFT)  and  the 
inverse  FFT  algorithms  for  permutohedral  aggregates  addressed  by  digital  indexing 
sets. 

Throughout,  assume  D,  X  €  Z+  fixed  and  let  n  =  2D+1  -  1.  Recall  that 
n  is  the  cardinality  of  a  Z)-dimensional,  A-level  permutohedral  aggregate.  We  will 
use  the  digital  indexing  sets  T\  and  T\  to  address  spatial  and  frequency  domain 
data,  respectively,  as  discussed  in  Section  5.4.  As  stated  earlier,  very  efficient  data 
processing  results  from  the  use  of  digital  indexing  sets  to  address  permutohedral 
aggregates  and  their  geometric  duals. 
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Define  the  map 

by 

(r,s)  =  exp[—27rf  (r,s)], 

where  Gx  is  a  set  of  coset  representatives  of  Ax  =  L/NXL,  G'x  is  a  set  of  coset 
representatives  of  A'x  =  L'X/(N*)XL'X,  and  (r,  s)  is  the  inner  product  of  r,  s  €  RD 
(see  Proposition  17);  (•,•)  is  a  bilinear  form  on  L  x  L'x.  For  a  6  L(A\),  define 
a  :  Gx  -»■  C  by  o(r)  =  «(r),  for  all  r  €  Gx;  for  a'  €  I(j4'a),  define  a'  :  GA  -+  C  by 
a'(s)  =  a'(s),  for  all  s  E  G'x.  If  g  €  T\  is  the  digital  address  of  r  €  GA,  and  /i  G  T'x 
is  the  digital  address  of  s  €  GA,  then  define  (g,  h)  =  (r,  5).  Hence,  in  the  context  of 
digital  addressing,  (•,  •)  is  a  map  Tx  x  T'x  — >  £/(nA). 

6.2     The  Discrete  Fourier  Transform  and  Inverse  DFT 

Let  a,fi  £  Z+  such  that  A  =  a  +  [3.    For  the  direct  transform,  define  ea  G 
Cixrj,ra  (a  row  vector  of  length  n2a)  by 

ea(v°qa)  =  (N'(f),V°)  6  U(nx),      q°  €  ra,t;Q  €  T'a. 

Let  e  =  ei.  Note  that  ex(hxgx)  =  (gx,hx)  and  e(hg)  =  (Nx-l(g),h).  For  the  inverse 
DFT,  define  e'a  e  CixrQr^  (a  row  vector  of  length  n2a)  by 


e'a(q°v")  =  (N0(q°),  v*)'1  =  (N'tf), »-)  €  £/(nA),      ga  €  TQ,  t,<*  €  I", 


Let  e'  =  e[.  Note  that  e'x(gxhx)  =  (gx,hx)  and  e'(gh)  =  (N^ig),^. 

Theorem  24  (Digital  DFT  for  Permutohedral  Aggregates)  IfD,\  €  Z+,n  = 
2D+1  -  1,  a  e  L(AA),  and  a'  =  jFA(a),  */*en 

.,  1 


a  = 


?} 


A/2 


(eA  0nA  of . 
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Proof.     The  discrete  Fourier  transform  of  a  €  L(A\)  is  determined  by 


nA/2 

r€GA 


Let  #A  €  Ta  be  the  digital  address  of  r  €  Ca  and  /iA  6  TA  the  digital  address  of 
s  €  G'x.  Then  equation  6.1  is  equivalent  to 


=   J72  E  ^(/*V  )<%A),    h"  e  r\. 

sA€rA 
In  terms  of  the  nA-product,  the  last  summation  can  be  written 


Q'  =  ^V2  (eA  ®"A  0)'  (6-2) 

to  yield  a  column  vector  of  length  nA.  □ 


The  nA-product  in  equation  6.2  requires  no  more  than  n2X  complex  multipli- 
cations. A  similar  proof  yields  the  corresponding  result  for  computation  of  the  digital 
inverse  DFT,  which  recovers  a  €  L(A\)  from  a'  €  L{A'X). 

Theorem  25  (Digital  Inverse  DFT  for  Permutohedral  Aggregates)  IfD,X  € 

Z+,n  =  2D+1  -  1,  a'  €  L(A'X),  and  a  =  T^l{a$)}  then 

a  =  ~^  (eA  ®«*  Q,y  ■ 

6.3     Coolev-Tukey  Factorization  Algorithms 

Let  a,p  €  Z+  such  that  A  =  a  +  (5.  For  the  Cooley-Tukey  (CT)  factorization 
of  the  direct  transform,  define  the  na  x  n0  matrix  Da<p  6  CpQ  XP   by 

D^iq^w1*)  =  (qa,wV)  e  U(nX),      qa  £  Ta,w0  e  V' 
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For  factorization  of  the  inverse  transform,  define  the  n0  x  n°  matrix  D'q  a  €  Cp'  xr0 
by 


VfiAv*,  <f)  =  (<T,  vfi)  €  U(nx),     w0  €  1%  gQ  6  TQ. 
Lemma  26  Ftt  L>,  A  €  Z+  and  let  a,/3  £  Z  suc/i  </iai  A  =  a  +  /?. 
/.  ///  6  T/j  and  q°  6  Ta,  tfien  Na{jpP)  +  qa  =  ^qa  €  TA; 
2.  ifva  <E  T'a  and  w0  €  T'p,  then  (N')0(va)  +  w0  =  vaw0  e  T'x. 
Proof.     Part  (1).  If  p0  =  p0_1  ■  ■  ■  plPo  and  q°  =  qQ.x  ■  ■  ■  qxqQ,  then 

pPff    =    P/3-i  ■  ■  ■  PiPoqa-i  ■  ■  •  <7i9o 

a+0-1  o-l 

=       Yl    Nk(Pk.a)  +  J2Nk(qk) 
k-a  k=0 

0-1  a-1 

k=0  k=0 

k=0  k=0 

Part  (2)  is  proved  similarly.  □ 

Lemma  27  Fix  Z),A  <E  Z+.  IfreL  and  s  €  L'x,  then  (Nx(r),s)  =  1. 

Proof.      We  have  Nx(r)  G  NXL.    Since  a  6  L'x  =  {NxL)',(Nx(r),s)  £  Z,  which 
implies  (NA(r),  5)  =  exp[-27ri(iVA(r), «)]  =  1.  □ 

Theorem  28  (Digital  CT  Factorization)  //  D,\,a,0  €  Z+  with  \  =  a  +  0, 
n  =  2D+1  -  1,  and  a  G  L(AX),  then 

e\  ©n»  a  =  ea  ©„o  (Z)aJj  o  (e^  ©n/S  a)) . 
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Proof.     The  n  -product  e\  ©n>  a  is  equivalent  to  the  summation 

£  eA(ftV)«(5A),     hx  e  r'A. 
3AerA 

If  gx  =  jPtf  e  TfiTa  =  TA  and  hx  =  vaw0  €  T'aY'0  =  PA,  then 

ca(AV)    =    (<A/>A) 

=     (jrqa,vawP) 

=    {NQ{p^)  +  9°,  (A^K)  +  «^)         (by  Lemma  26) 
.     =    (^(/),(Af*)/JK))-(iVQ(p/3),^)-(9a,(^)>Q)).(^,^) 

=  (wV),o  •  (^a(/)X)  •  (^V),»a)  •  (<rX) 

the  last  equality  by  Lemma  27.  Hence, 

£  ex(hY)atf)    =     EE  eAKu;VV)a(pV)  (6.3) 

5AerA  7QerQp/56r/3 

=     £  C(»Y)  [^(fy^fl-y)]  ,         (6.4) 
9Qera 

where 

CM*,*?)*    £   e/»(«V)«G»V).  (6.5) 

Summation  6.4  is  equivalent  to  eQ  ©na  (jDajS  o  Cp)  and  equation  6.5  is  equivalent  to 
Cp  =  ep  ©„^  a.  Hence, 

ex  0„a  a    =     ea  ©„<,  (Dat0  o  C» 

=     ea  0„a  (Z)Q>/8  o  (e/j  ©„*>  a)), 
a  row  vector  of  length  nx.  rj 

The  inside  ^-product  in  the  CT  factorization  algorithm  ea©na  {Da<po{ep®np 
a))  requires  no  more  than  n2p  complex  multiplications;  the  Hadamard  product  takes 
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no  more  than  nA  complex  multiplications.  Hence,  the  total  number  of  complex  mul- 
tiplications required  by  the  CT  factorization  algorithm  is  bounded  by  n2a  +  nw  +  nx. 

Remark  6   Theorems  24  and  28  imply 

a'  =  —rr;  {ea  8*«  (Da,p  o  (ep  ®nt>  a)))' . 

Arguments  similar  to  those  in  the  proof  of  Theorem  28  establish  the  digital 
inverse  CT  factorization  algorithm. 

Theorem  29  (Digital  Inverse  CT  Factorization)  If  D,\,a,p  e  Z+  with  A  a 
a-r{3,n  =  2D+1  -  1,  a'  e  L(A'X),  and  a  =  7^{a'),  then 

e\  ©„a  a'  =  e'0  ®n,  (D'0a  o  (e'a  ©„*  a'))  . 

Remark  7  Theorems  25  and  29  imply 

■■J^W  ®n?  (D'0>a  o  (e'a  ©na  a')))1 . 

Theorems  28  and  29  are  digital  versions  of  the  CT  factorization  for  permutohe- 
dral  aggregates.  However,  the  next  two  lemmas,  which  use  variations  of  the  matrices 
Datf3  and  Dp  ,  respectively,  yield  digital  FFT  algorithms  (direct  and  inverse)  for 
permutohedral  aggregates.  Define  Ta+itp-i  €  Cra+lXr_    by 

Tt+ij-t(pu4*,w0-1)  =  {Na{PM\wP-1)  6  U(nx),     pMqa  6  r^t V1  6  ^_, 

and  define  I^H,*-!  €  Cr^+1  xr,-,  by 

In  the  signal  processing  literature,  the  entries  in  matrices  Da<p,  D'p  a,Ta+itp-1,  and 
Tp+la_l  are  known  as  twiddle  factors. 
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Lemma  30  LetD,\  €  Z+,n  =  2D+1-1,  and  a  €  L(AX).  If  P  €  Z+  with  2  <  /?  <  A, 

e^3  ©n^  «  =  e  0n  (Ta+ijj-i  o  (e^_i  ©n/5-i  a)) , 
where  a  =  A  —  /?. 

Proof.    The  n^-product  e/?  ©„?  a  is  equivalent  to  the  summation 

£  *C*V)t(pV).    ?aera,^er;. 

p"er5 
If/  =  /-xpm  €  r^-,r,  =  T^  and  u/  =  wMw0-1  €  rjl^.,  =  r^,  then 

^(wV)  =   (ivV)V) 

=  (^(/-'pm),^-1) 

=  (/va+ V-1)  +  n"(pm),  (tr-f-^WM)  +  w0-1) 

=    (Na+1(pP-1)1(N'f-1(wM))-(Na+1(^-1)t^-1) 

=     1  •  e^u^-y-1)  •  e(u;A/pM)  •  T^^^cf  V-1), 
the  last  equality  by  Lemma  27.  Hence, 

5]  e/?(iA>y<n  =  E    E  ^AX-y-v)ay-v<n 

pmgi\ 

•e(^A/PA/),  (6.6) 
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where 

C0-1(pMqa,We-1)=       Y,      ^(t^-y-'MAW*)-  (6-7) 

Summation  6.6  is  equivalent  to  e  ©„  (To+i^-i  o  Cp-\)  and  equation  6.7  is  equivalent 
to  Cp-i  =  ep_i  ©„/s-i  a,  a  matrix  of  dimensions  n°  x  n^.  Hence, 

=    e  0„  (Tb+i^-i  o  (e/j-i  ©n/i-i  a)). 


Similar  arguments  establish  the  corresponding  inverse  digital  algorithm. 

Lemma  31  LetD,X  €  Z+,n  =  2D+1-1,  and  a'  6  L(A'A).  //a  €  Z+  with  2  <  a  <  X, 
then 

e'a  0„-  a'  =  e'  ©n  (2J+i*-i  o  «_,  0n«_,  a'))  , 
where  (3  —  X  —  a. 


6.4     The  Fast  Fourier  Transform  and  Inverse  FFT  Algorithms 

The  following  algorithm  is  a  radix-n,  decimation-in-space,  fast  Fourier  trans- 
form for  permutohedral  aggregates  addressed  by  digital  indices. 

Theorem  32  (Digital  FFT  for  Permutohedral  Aggregates)  LetD,X  6  Z+,n 
2D+1-l,o€  L{AX),  and  a'  =  Fx{a).  If 

Cx  =  e  0n  a,  (6.8) 

and 

C0  =  e®n(Tx-0+i,0-ioC0-i),  for  (3  =  2,3,. ..,X,  (6.9) 
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then 


1 

a'  = 


>  -     *    cl 


Proof.  We  claim  Cp  =  tp  ©n/s  a,  for  1  <  /3  <  A.  The  statement  is  true  for  /?  =  1  by 
hypothesis.  Arguing  by  mathematical  induction,  if  C7  =  e7  ©n7  a,  for  1  <  7  <  A  —  1, 
then 

C7+1     =    e  ©n  (7a_7)7  o  C7) 

=     e  ©n  (Ta_7i7  o  (e7  ©n7  a)) 

=     e7+i  ©„7+i  a,  (by  Lemma  30) 

for  1  <  7  <  A  -  1.  Hence,  Cp  =  ep  ©n,s  a,  for  1  <  /?  <  A.  In  particular,  C\  =  e>  ©„*  a, 
which  implies  a'  =  n~x/2  ■  C{,  by  Theorem  24.  D 

A  similar  argument  proves  the  corresponding  digital  inverse  FFT  for  permu- 
tohedral  aggregates,  a  radix-n,  decimation-in-frequency  algorithm. 

Theorem  33  (Digital  Inverse  FFT  for  Permutohedral  Aggregates)  Let  D,  A  6 

Z+,n  =  2D+1  -l,o'€  L{A\),  and  a  =  F^{a').  If 

C/j  =  e  ©n  a  , 

and 

^  =  e'©„(n_Q+1,0_1oC:_1),/ora  =  2,3,...,A, 


n 


A/2' 


ice 


Computation  of  the  DFT  of  a  function  a  £  L{AX)  consisting  of  nx  latt 
points  or  cells  requires  A  stages,  as  specified  by  equations  6.8  and  6.9.  Computation 
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of  the  n-product  in  each  of  these  equations  requires  no  more  than  n  +1  complex 
multiplications.  The  Hadamard  product  in  each  of  the  A  —  1  stages  specified  by 
equation  6.9  requires  at  most  nA  complex  multiplications.  Hence,  the  total  number  of 
complex  multiplications  required  to  compute  the  DFT  of  a  using  the  FFT  algorithm  is 
bounded  by  AnA+1  +  (A  —  l)nA.  The  algorithm,  therefore,  has  complexity  0(M  log  M), 
where  M  —  nx.  Furthermore,  since  the  FFT  algorithm  (Theorem  32)  is  written 
entirely  in  terms  of  the  generalized  matrix  product  and  Hadamard  multiplication,  no 
explicit  reordering  of  the  data  is  required.  Data  reordering  is  a  common  feature  of 
standard  FFT  algorithms.  Similar  remarks  apply  to  the  computation  of  the  digital 
inverse  FFT  algorithm  (Theorem  33). 


CHAPTER  7 
SUMMARY  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


In  this  dissertation  we  showed  how  the  additive  structure  of  permutohedral 
aggregates  can  be  exploited  to  derive  fast  algorithms  for  the  computation  of  the 
discrete  Fourier  transform  of  functions  whose  domains  are  aggregates  of  any  given 
level  in  any  given  dimension.  In  these  algorithms,  the  data  is  addressed  in  terms  of 
an  efficient  digital  addressing  scheme  known  as  generalized  balanced  ternary  (GBT) 
addressing.  The  algorithms  described  here  are  not  based  on  the  Chinese  remainder 
theorem,  which  is  used  to  derive  the  fast  Fourier  transform  (FFT)  for  functions  de- 
fined on  rectangular  arrays.  Furthermore,  we  defined  the  generalized  matrix  product 
for  permutohedral  aggregates  and  stated  the  FFT  and  the  inverse  FFT  for  permuto- 
hedral aggregates  entirely  in  terms  of  digital  addressing  and  the  generalized  matrix 
product.  The  notation  described  in  the  text  renders  very  succint  algorithms  that  are 
amenable  to  computer  implementation. 

Note,  however,  that  we  have  not  directly  exploited  the  multiplicative  structure 
of  permutohedral  aggregates.  A  number  of  algorithms  exist  that  compute  the  discrete 
Fourier  transform  for  domains  that  possess  a  multiplicative  structure  (see  Tolimieri, 
An,  and  Lu  [20]).  Hence,  we  propose  the  following  problems  for  further  research: 

1 .  characterization  of  ideals  and  quotients  of  permutohedral  aggregates  in  terms 
of  GBT  addressing; 

2.  derivation  of  the  discrete  Fourier  transform  of  periodic  and  decimated  functions 
defined  on  permutohedral  aggregates; 
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3.  derivation  of  DFT  algorithms  that  exploit  the  multiplicative  structure  of  per- 
mutohedral  aggregates; 

4.  extension  of  the  definition  of  the  generalized  matrix  product  for  permutohedral 
aggregates  to  more  general  addressing  schemes; 

5.  description  of  self-similar  tiles  derived  from  permutohedral  aggregates  and  deriva- 
tion of  a  discrete  wavelet  transform  and  a  multiresolution  analysis  for  such  tiles. 


REFERENCES 


[1]  M.  An,  I.  Gertner,  M.  Rofheart,  and  R.  Tolimieri,  Discrete  fast  Fourier  transform 
algorithms:  a  tutorial  survey,  Advances  in  Electronics  and  Electron  Physics,  Vol. 
80,  Academic  Press,  Inc.,  New  York,  NY,  1991,  pp.  1-67. 

[2]  L.  Auslander  and  R.  Tolimieri,  Is  computing  with  the  Fourier  transform  pure 
or  applied  mathematics?,  Bulletin  (New  Series)  of  the  American  Mathematical 
Society,  Vol.  1,  No.  6,  November  1979,  pp.  847-897. 

[3]  C.  Bandt,  Self-similar  sets  5.  Integer  matrices  and  fractal  tilings  of  Rn,  Pro- 
ceedings of  the  American  Mathematical  Society,  Vol.  112,  No.  2,  June  1991,  pp. 
549-562. 

[4]  J.  H.  Conway  and  N.  J.  A.  Sloane,  Sphere  packings,  lattices,  and  groups, 
Springer- Verlag,  New  York,  NY,  1988. 

[5]  J.  W.  Cooley  and  J.  W.  Tukey,  An  algorithm  for  the  machine  calculation  of 
complex  Fourier  series,  Math.  Comp.,  Vol.  19,  No.  2,  1965,  pp.  297-301. 

[6]  C.  L.  DeVito,  Functional  analysis  and  linear  operator  theory,  Addison- Wesley 
Publishing  Company,  New  York,  NY,  1990. 

[7]  D.  E.  Dudgeon  and  R.  M.  Mersereau,  Multidimensional  digital  signal  processing, 
Prentice  Hall,  Englewood  Cliffs,  NJ,  1984. 

[8]  L.  Gibson  and  D.  Lucas,  Spatial  data  processing  using  generalized  balanced 
ternary,  Proceedings  of  the  IEEE  Conference  on  Computer  Vision  and  Pat- 
tern Recognition,  June  1982,  pp.  566-571. 

[9]  R.  C.  Gonzalez  and  R.  E.  Woods,  Digital  image  processing,  Addison-Wesley, 
Reading,  MA,  1992. 

[10]  T.  W.  Hungerford,  Algebra,  Springer- Verlag,  Berlin,  1980. 

[11]  M.  T.  Heideman,  D.  H.  Johnson,  and  C.  S.  Burrus,  Gauss  and  the  history  of  the 
fast  Fourier  transform,  IEEE  ASSP  Magazine,  October  1984. 

[12]  J.  Kappraf,  Connections:  the  geometric  bridge  between  art  and  science,  McGraw- 
Hill,  Inc.,  New  York,  NY,  1991. 

[13]  W.  Z.  Kitto,  A.  Vince,  and  D.  C.  Wilson,  An  isomorphism  theorem  between  the 
p-adic  integers  and  a  ring  associated  with  a  tiling  of  N -space  by  permutohedra, 
Discrete  Applied  Mathematics  52,  1994,  pp.  39-51. 


69 


70 


14]  D.  Lucas,  A  multiplication  in  N -space,  Proceedings  of  the  American  Mathemat- 
ical Society  74,  1979,  pp.  1-8. 

15]  R.  M.  Mersereau  and  T.  S.  Speake,  The  processing  of  periodically  sampled  mul- 
tidimensional signals,  IEEE  Transactions  on  Acoustics,  Speech,  and  Signal  Pro- 
cessing, ASSP-31,  No.  1,  1983,  pp.  188-194. 

16]  D.  P.  Petersen  and  D.  Middleton,  Sampling  and  reconstruction  of  wave-number 
limited  functions  in  n-dimensional  euclidean  spaces,  Information  and  Control, 
Vol.  5,  1962,  pp.  279-323. 

17]  G.  X.  Ritter,  Heterogeneous  matrix  products,  in  Image  algebra  and  morpholog- 
ical image  processing  II,  Vol.  1568  of  Proceedings  of  SPIE,  San  Diego,  CA,  July 
1991,  pp.  92-100. 

18]  G.  X.  Ritter  and  J.  N.  Wilson,  Handbook  of  computer  vision  algorithms  in  image 
algebra,  CRC  Press,  Boca  Raton,  1996. 

19]  G.  Strang,  Wavelet  transforms  versus  Fourier  transforms,  Bulletin  of  the  Amer- 
ican Mathematical  Society,  Vol.  28,  No.  2,  April  1993,  pp.  288-305. 

20]  R.  Tolimieri,  M.  An.,  and  C.  Lu,  Algorithms  for  discrete  Fourier  transform  and 
convolution,  Springer- Verlag,  New  York,  1989. 

21]  R.  Tolimieri,  M.  An.,  and  C.  Lu,  Mathematics  of  multidimensional  Fourier 
transform  algorithms,  2nd  ed.,  Springer- Verlag,  New  York,  1997. 

22]  A.  Vince,  Replicating  tesselations,  SIAM  Journal  of  Discrete  Mathematics,  Vol. 
6,  No.  3,  August  1993,  pp.  501-521. 

23]  A.  B.  Watson  and  A.  J.  Ahumada,  Jr.,  A  hexagonal  orthogonal-oriented  pyramid 
as  a  model  of  image  representation  in  the  visual  cortex,  IEEE  Transactions  on 
Biomedical  Engineering,  Vol.  36,  No.  1,  1989,  pp.  97-106. 

24]  H.  Zhu,  The  generalized  matrix  product  and  its  applications  in  signal  processing, 
Ph.D.  Thesis,  1993. 

25]  H.  Zhu  and  G.  X.  Ritter,  The  generalized  matrix  product  and  the  wavelet  trans- 
form, Journal  of  Mathematical  Imaging  and  Vision,  Vol.  3,  No.  1,  1993,  pp. 
95-104. 

26]  H.  Zhu  and  G.  X.  Ritter,  The  p-product  and  its  applications  in  signal  processing, 
SIAM  Journal  of  Matrix  Analysis  and  Applications,  Vol.  16,  No.  2,  1995,  pp. 
579-601. 


BIOGRAPHICAL  SKETCH 

Jaime  L.  Zapata  was  born  in  Bogota,  Colombia.  He  received  the  Bachelor 
of  Science  degree  in  engineering  with  High  Honors  from  the  University  of  Florida  in 
August  1985,  and  the  Master  of  Science  degree  in  biomedical  engineering  from  the 
University  of  Miami  in  August  1986.  From  1986  to  1990,  he  worked  as  a  software 
engineer  with  Coulter  Electronics,  a  major  biomedical  engineering  firm  in  Hialeah, 
Florida.  His  research  interests  include  image  algebra,  image  processing,  wavelets  and 
multiresolution  analysis,  and  bioinformatics  and  biocomputing. 


71 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  itxonforms  to 
acceptable  standards  of  scholarly  presentation  and  is  fully  adequate^  ii)  scope  and 
quality,  as  a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


rd  X.  Ritter,  Chairman 
of  Mathematics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it  conforms  to 
acceptable  standards  of  scholarly  presentation  and  is  fully  adequate,  in  scope  and 
quality,  as  a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


David  C.  Wilson 
Professor  of  Mathematics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it  conforms  to 
acceptable  standards  of  scholarly  presentation  and  is  fully  adequate,  in  scope  and 
quality,  as  a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


\iwt Ju    If 


A-*L 


Murali  Rao 

Professor  of  Mathematics 


I  certify  that  I  have  read  this  study  and  that  in  my  opinion 
acceptable  standards  of  scholarly  presentation  and  is  fully  adequ 
quality,  as  a  dissertation  for  the  degree  of^o^ctor  of  Philosophy. 


£ 


it  conforms  to 
te,  in  scope  and 


JameSi  E.  Keesling 
/Professor  of  Mathematics 

I  certify  that  I  have  read  this  study  and  that  in  my  opinion  it  conforms  to 
acceptable  standards  of  scholarly  presentation  and  is  fully  adequate,  in  scope  and 
quality,  as  a  dissertation  for  the  degree  of  Doctor  of  Philosophy. 


Jo>ieph  N.  Wilson 

Assistant  Professor  of  Computer  and 
Information  Science  and  Engineering 


This  dissertation  was  submitted  to  the  Graduate  Faculty  of  the  Department  of 
Mathematics  in  the  College  of  Liberal  Arts  and  Sciences  and  to  the  Graduate  School 
and  was  accepted  as  partial  fulfillment  of  the  requirements  for  the  degree  of  Doctor 
of  Philosophy. 

May  2000 


Dean,  Graduate  School 


UNIVERSITY  OF  FLORIDA 

illlllllllllllllllllllll 

3  1262  08555  3617 


